NASA Technical Memorandum 100509 


POSTBUCKLING OF LAMINATED ANISOTROPIC PANELS 


(HASA-TH-10G5C9) POSlfiDCKIlliC Cf 1AMIKA1ED N88-16G11 

i AISC1ECEIC PAHELE (KASA) 182 p CSCL 20K 

Unclas 

G3/39 0103611 


Glenda L. Jeffrey 


October 1987 


NASA 

National Aeronautics and 
Space Administration 

Langley Research Center 

Hampton Virginia 23665 


Acknowledgements 


Many people at Langley Research Center have been indispensable for their assistance 
in completing the present work. The author expresses her gratitude first to Dr. Ahmed 
K. Noor, without whom this work would not have been possible. She would also like to 
thank the employees of the Structural Mechanics Branch for time spent in many useful 
discussions. Special thanks is extended to Dr. James H. Starnes, who helped to guide this 
work, and to Dr. Norman F. Knight, Jr., whose help in running STAGS was invaluable. 
George Johnson and Jim Kiss, the technicians involved in the experimental work, were 
excellent; the author could not have asked for better. Gratitude is extended to Jeanne M. 
Peters for her programming assistance, and to Carl M. Andersen for many useful discus- 
sions. The author is also thankful to her supervisors, Brantley R. Hanks and Dr. Jerrold 
M. Housner, for their tolerance during her first months of employment. Finally, the au- 
thor is especially grateful for the knowledge, tolerance, moral support, and computational 
resources volunteered by the personnel of the Computational Structural Mechanics Group. 


PRECEDING PAGE BLANK NOT FILMED 


m 


Table of Contents 


Abstract i 

Acknowledgements iii 

Table of Contents v 

List of Tables ix 

List of Figures xi 

List of Symbols xv 

1. Introduction 1 

1.1. Overview 1 

1.2. Review of Previous Pertinent Work 3 

1.2.1. Global Approximation Function Techniques 3 

1.2.2. Reduced Basis Work 4 

1.2.3. Experimental Work 5 

1.3. Objectives and Scope 6 

2. Problems Involving Bending-Extensional Coupling and Anisotropy . S 

2.1. Mechanical Couplings S 

2.2. Symmetries Exhibited by Anisotropic Panels 13 

2.3. Governing Differential Equations of an Anisotropic Plate 

And Bending-Extensional Coupling 16 

3. Development of the Modified Rayleigh-Ritz Technique 18 

3.1. Review of the Rayleigh-Ritz Technique 18 

3.1.1. Development of the Exact Energy Expression IS 

3.1.2. Selection and Use of Approximation Functions 20 

3.1.3. Generation of the Nonlinear Equations in Xj And the Recursion 

Formulas for the Newton-Raphson Iterative Procedure 23 

3.1.4. Symmetrization of the Nonlinear Arrays 24 

3.1.5. Generation of Load Vector P 24 

3.2. Enforcement of Boundary Conditions Using 

A Penalty Function Approach 25 

3.3. Operator Splitting 27 

v 


PAGE — Tp?* INTtNTKMU-y BLANK 


.^/PRECEDING PAGE BLANK NOT FILMED 



‘i ..‘{.I Basic Concepts 28 

.'>.'>.2. A 1 1 ;iy Symmetri/.ation When lhi.rl.il, ioning is Used 32 

3.4. Application of Reduced Basis Techniques in Conjunction 

With Operator Splitting 34 

3.4.1. Generation of Basis Vectors 35 

3.4.1. a. Linear Problem 35 

3.4.1. b. Nonlinear Problem 37 

3.4.2. Generation and Solution of the Reduced Equation System 41 

3.4.2. a. Linear Problem 41 

3.4.2. b. Nonlinear Problem 41 

4. Numerical Studies of Modified Rayleigh- Ritz Technique 44 

4.1. Comments on Problem Selection 44 

4.2. Linear Results 44 

4.2.1. Transversely Loaded Orthotropic Panel 47 

4.2.2. Transversely Loaded Anisotropic Panel With 

Bending- Extensional Coupling 49 

4.2.3. Axially Loaded Orthotropic Panel 50 

4.2.4. Axially Loaded Anisotropic Panel With 

Bending-Extensional Coupling 52 

4.3. Effect of the Penalty Parameter on Accuracy of Linear Results 53 

4.4. Linear Results Using Reduced Basis Techniques 55 

4.4.1. Accuracy in Approximating the Full System 57 

4.4.2. Grammian Results 58 

4.5. Nonlinear Results 59 

4.5.1. Orthotropic Panel 63 

4.5.2. [±30/90] 5 Panel 64 

4.5.3. [±50/ 35] 5 Panel 65 

5. Experimental and Numerical Studies of 

Buckling and Postbuckling Behavior of Panels 67 

5.1. Characteristics of Test Specimens 67 

5.2. Description of the Experiment 68 

5.2.1. Test Setup, Measurement Techniques, and Procedure 68 

5.2.2. Observations Regarding Experimental Techniques 73 

5.3. Finite Element Model 74 

5.3.1. Discretization and Boundary Conditions 74 

5.3.2. Solution Strategy 76 

5.4. Comparison of Experimental and Finite Element Results 77 

5.4.1. Presentation and Discussion of Buckling Load 

And End Shortening Results 77 

5.4.1. a. Buckling Load Results and Discussion 78 

5.4.1. b. Postbuckling Results and Discussion 81 

5.4.2. Presentation and Discussion of Out-of-Plane Displacement Results . 84 

vi 



5.4.3. Presentation and Discussion of Strain Results 86 

5.5. Failure Initiation and Development 89 

5.5.1. Qualitative Results 89 

5.5.2. Quantitative Results 91 

5.6. Summary of Observations 91 

6. Observations and Recommendations 97 

6.1. Observations 97 

6.1.1. Observations Related to the Modified Rayleigh- Ritz Technique ... 97 

6. 1.1. a. Observations related to linear results 97 

G.l.l.b. Observations related to nonlinear results 98 

6.1.2. Observations Based on Experiments 99 

6.1.3. Observations Based on Finite Element Modeling Efforts 100 

6.2. Recommendations 102 

6.2.1. Recommendations Related to the Modified Rayleigh- Ritz Technique 102 

6.2.2. Recommendations Related to Experimental Work 102 

References 104 

Appendices 

A. Array Symmetrization Using Only Quadratic Energy Terms . . 106 

B. Experimental Moire Fringe Photographs 108 

C. C-Scans of Experimental Panels 116 

D. STAGSC-1 Runstreams 123 

Figures 130 



List of Tables 


2.1. Orthotropic Panel Symmetry States 15 

4.1. Rayleigh- Ritz Analysis Functions 46 

4.2. Linear Results for [O 4 / 9 O 43 S, Transversely Loaded Panel 48 

4.3. Linear Results for [±50/35] s , Transversely Loaded Panel 50 

4.4. Linear Results for [04/904]s, Axially Loaded Panel 51 

4.5. Linear Results for [±50/35] 5 , Axially Loaded Panel 52 

4.6. Penalty Results - Overall Penalty Parameter 54 

4.7. Variation of Results with Number of Basis Vectors Used 57 

4.8. Variation of the Grammian with Number of Basis Vectors Used 58 

5.1. Material Properties for Hercules AS4-3502 Graphite/Epoxy Tape 68 

5.2. Constitutive Matrices for [±50/35]s and [±30/90]5 Panels, Thickness = 0.08 in. 69 

5.3. Measurements of Tested Panels 70 

5.4. Finite Element Model Discretization 74 

5.5. Grid Convergence, [±30/90]s 4-inch-wide Panel 75 

5.6. Buckling Loads, Pounds 79 

5.7. Linear Bifurcation Buckling Results, Pounds 80 

5.8. Slope of Initial Postbuckling Curve 81 

5.9. Postbuckling Stiffness Degradation 82 

5.10. Transverse Displacement at Buckling 85 

5.11. Summary of Twist and Symmetry Results 86 

5.12. Load and End Shortening at Failure 92 

dRECEDiNG PAGE BLANK NOT FILMED 

ix 

FK£ - Jp T WTtHHONAHY BUNK 


List of Figures 


B.l. Moire Fringe Photographs for [±30/90] 5 4 Inch-Wide Panel 109 

B.2. Moire Fringe Photographs for [±30/90]s 5.5 Inch-Wide Panel 110 

B.3. Moire Fringe Photographs for [±30/90] 5 7 Inch- Wide Panel Ill 

B.4. Moire Fringe Photographs for [±50/35] 5 4 Inch- Wide Panel 112 

B.5. Moire Fringe Photographs Illustrating Mode Change 

Exhibited by [±50/35] 5 4 Inch- Wide Panel 113 

B.6. Moire Fringe Photographs for [±50/35]s 5.5 Inch- Wide Panel 114 

B. 7. Moire Fringe Photographs for [±50/35]s 7 Inch-Wide Panel 115 

C. l. C-Scan of [±30/90] 5 4 Inch- Wide Panel 117 

C.2. C-Scan of [±30/90]s 5.5 Inch-Wide Panel 118 

C.3. C-Scan of [±30/90] 5 7 Inch- Wide Panel 119 

C.4. C-Scan of [±50/35] 5 4 Inch- Wide Panel 120 

C.5. C-Scan of [±50/35]s 5.5 Inch- Wide Panel 121 

C.6. C-Scan of [±50/35] s 7 Inch- Wide Panel 122 

2.1. Conventions Used in Theoretical Development 130 

2.2. Effect of Shear-Extensional Coupling on Panel Subject To 

Uniform Axial Loading 131 

2.3. Effect of Shear-Extensional Coupling on Postbuckled Shape 131 

2.4. Method for Observing Twisting Behavior During an Experiment 132 

4.1. Schematics for Linear Problems 133 

4.2. Schematic of Mixed Formulation Finite Element Models 134 

4.3. Linear Results for Simply Supported, Transversely Loaded [04/904] s Panel . . 135 

4.4. Linear Results for Simply Supported, Transversely Loaded [±50/35] 5 Panel . 136 

4.5. Linear Results for Clamped/Simply Supported, Axially Loaded [04/904] s Panel 137 


4.6. Linear Results for Clamped/Simply Supported, Axially Loaded [±50/35]5 Panel 138 


xi 


PRECEDING PAG S 8LA&X NOT FILMED 


KU law ALLY BUlQi 


4.7. Contour Plots of Linear Solution for Axially Loaded 

Clamped/Simply Supported [±50/35]5 Panel 139 

4.8. End Shortening Behavior, [ 04 / 90 4 ] s Panel 140 

4.9. Center Transverse Displacement Behavior, [0 4 /90 4 ]s Panel 141 

4.10. Normalized Contour Plots of Transverse Displacement, [0 4 /90 4 ] s Panel . . . 142 

4.11. End Shortening Behavior, [±30/90] 5 Panel 143 

4.12. Center Transverse Displacement Behavior, [±30/90]5 Panel 144 

4.13. Normalized Contour Plots of Transverse Displacement, [±30/90]s Panel . . . 145 

4.14. End Shortening Behavior, [±50/35]s Panel 146 

4.15. Center Transverse Displacement Behavior, [±50/35 ] 5 Panel 147 

4.16. Normalized Contour Plots of Transverse Displacement, [±50/35]s Panel . . . 148 

5.1. Experiment Schematic and Photograph 149 

5.2. Support Fixture Schematic 150 

5.3. Bending as Indicated by Back-to-back Strain Gage Data 151 

5.4. Sample C-Scan 152 

5.5. Comparison of Experimental Setup Technique for Isotropic 

And Anisotropic Panels 153 

5.6. Undeformed Geometry Plot of STAGSC-1 Finite Element Model 

for 5.5-inch-wide [±50/35]s Panel 154 

5.7. Typical Panel End-shortening Plot 155 

5.8. End Shortening Results 156 

5.9. Estimation of Stiffness Degradation 157 

5.10. Near-Maximum Transverse Displacement Results 158 

5.11. Experimental Warping Results 159 

5.12. STAGSC-1 Warping Results 160 

5.13. Experimental Symmetry Results 161 

5.14. STAGSC -1 Symmetry Results 162 

xii 



5.15. Representative Experimental Surface Strain Results 1 03 

5.16. Progression of Transverse Displacement Pattern, [±50/35]s, 7-inch- wide panel 164 

5.17. Progression of Transverse Displacement Pattern, [±30/90] 5 , 7-inch-wide panel 165 


5.18. Progression of Transverse Displacement Pattern, All Panels 166 

5.19. Failure Locations 167 

5.20. Strains Across 5.5-inch-wide Panels at Load Near Failure 168 

5.21. Photograph of Shear Failure, [±50/35]5 7-inch-wide panel 169 


xm 


List of Symbols 


a 

b 

A 

B 

C 

c 

D 


F 

F 

F 

G 

G 

h 

X 

K 

K 

K 

M 

M 

m 

M 

N 


Longer panel dimension 

Shorter panel dimension 

Membrane stiffness coefficients 

Bending-extensional coupling stiffness coefficients 

Constitutive matrix 

Transverse shear stiffness coefficients 

Bending stiffness coefficients 

Error measure used in Newton-Raphson iterations 

Vector of approximation functions representing displacement i 

and containing j terms 

Vector containing all Fj for all displacements 
Symmetrized quadratic stiffness matrix which is cubic in Fj 
Unsymmetrized quadratic stiffness matrix which is cubic in Fj 
Symmetrized cubic stiffness matrix which is quartic in Fj 
Unsymmetrized cubic stiffness matrix which is quartic in Fj 
Constraint function; also panel thickness (clear from context) 
Total number of even state approximation functions 
Symmetrized linear stiffness matrix which is quadratic in Fj 
Unsymmetrized linear stiffness matrix which is quadratic in Fj 
Reduced linear stiffness matrix 
Moment resultant 

Total number of odd state approximation functions 
Total number of approximation functions 
Gram matrix 
Stress resultant 


xv 

PRECEDING PAGE BLANK NOT FILMED 


.WTEtf7rON.Ai.lY etAN* 



n 

M 


P 

P 


p 

Q 

r l 


r nl 


R 

r 

U 


Ui 


Ui 


«2 e 

^ext 

w 

W c 

Xi 

x 'i 

X 

X 

6 

8ij 


Total number of degrees of freedom in nonlinear equation system 

Number of basis vectors to be used 

Panel buckling load 

Load vector 

Reduced load vector 

Load magnitude 

Shear stress resultants 

Matrix containing derivatives of the Fj which constitute 

linear contributions to strain displacement relations 

Matrix containing derivatives of the Fj which constitute 

nonlinear contributions to strain displacement relations 

Residual vector 

Penalty parameter 

Total energy of a structure 

Strain energy terms involving generalized displacements raised to power i 

Displacement in direction i 

Value of i/2 at midpoint of long panel edge 

Work of external forces 

Transverse displacement 

Center transverse displacement 

Spatial coordinate i 

Vector of generalized coordinates representing displacement i 
and containing j terms 

Vector containing all Xj for all displacements 

Vector amplitudes associated with basis vectors contained in T 

End shortening 

Kronecker delta, equal to one when i = j and zero otherwise 


xvi 



8U 


£ 



V 

K 

A 
4>i 
P 2 

a 

6 

r 

n 


Variation of U 
Strain 

Linear contribution to strain-displacement relations 
Nonlinear contribution to strain-displacement relations 
Path parameter 
Curvature 

Tracing constant associated with anisotropic stiffness terms 
Rotation about axis £2 
Rotation about axis xi 
Stress 

Angle principle material axis of orthotropic lamina makes with x\ axis 
Matrix whose columns are basis vectors 
Penalty function 


Subscripts 

0 Evaluation at A = 0 

a, b 1, . . . , 8 

a, ft Dummy indices 

c, d 1, . . . , Af 

l, J, K, L 1, . . . ,m 

b j 1 ^ > l 1, . . . , I 

m, n, o, p 1 , . . . , M. 

P, Q, R, S, T Range is over number of functions approximating 

« 2 i <f>u P 21 and w, respectively 

All repeated indices denote summation unless otherwise noted. 


XVII 



o 


Orthotropir material stiffness coefficients 

Anisotropic (nonorthotropic) material stiffness coefficients 

Even and odd state functions, respectively 


A 

e, o 



XVlll 



Chapter 1 
Introduction 

1.1 Overview 

The potential of aeroelastic tailoring methods involving laminated filamentary com- 
posite aircraft wing skins has recently been exploited for two primary purposes: 

1. Achievement of performance goals by causing changes in wing shape for specified 

aerodynamic loads. 

2. Elimination of instabilities, such as wing divergence, from otherwise desirable designs. 
Most tailoring efforts have centered around obtaining desired twisting and bending be- 
havior through careful choice of laminate stacking sequences. Such tailoring is possible 
because many different mechanical couplings can be induced in a laminate through the use 
of certain well-defined stacking sequence categories. 

Two recent and frequently cited design applications employing aeroelastic tailoring 
have been the HiMAT remotely piloted vehicle [1] and X-29 fighter plane design efforts. The 
HiMAT (Highly Maneuverable Aircraft Technology) aircraft, a 0.5 scale research vehicle, 
is an example of the use of aeroelastic tailoring to achieve a performance goal. In this 
case, design constraints required that the aircraft’s wings attain a certain twist in response 
to certain operational loads. The X-29, on the other hand, is an example of the use of 
aeroelastic tailoring to eliminate instabilities from an otherwise feasible design. The wings 
of the X-29 were tailored in such a way that the divergence inherent in the forward-swept 
wing was avoided [2]. 

It is important to note that, even though bending/ twisting coupling is the interaction 
usually desired for aeroelastic tailoring applications, the [±50/35] 5 laminate used in the 
HiMAT aircraft exhibited all possible types of mechanical couplings. (See Chapter 2 for a 
discussion of all the different couplings.) Apparently, the analysis only included the effects 
of bending- and extension-twisting coupling, which might have accounted for some of the 
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test /analysis discrepancies discussed in [ 1 ]. Such an omission is perhaps understandable, 
in light of the dearth of information regarding the analysis of laminates simultaneously 
possessing all the couplings present in the [±50/35 ] 5 laminate. Almost all the literature 
found by the author deals with laminates in which one or more couplings are not present. 
It is therefore desirable to compile theoretical and experimental results dealing with such 
laminates. 

A second area which has received much attention is the practical use of laminated fil- 
amentary composites in the postbuckling range. Although metallic aircraft substructures 
having buckling loads within the operational load range of the aircraft have been in use 
for some time, composite aircraft structures have typically been designed by considering 
buckling to be synonymous with failure [3]. It has only been fairly recently that investiga- 
tors have shown that buckling need not be such a limitation, since considerable stiffness is 
often available at loads equal to many times the buckling load of a composite structure. 
(See, for example, [4] and [3].) Efficient use of this additional stiffness makes possible the 
design and construction of structures which are significantly lighter than metal structures 
of equivalent strength. 

The possible future use of lighter, buckled designs in aeroelastically tailored substruc- 
tures obviates the need for more complete information regarding both the pre- and post- 
buckling behavior of anisotropic composite laminates with bending-extensional coupling. 
Although considerable literature exists regarding postbuckling behavior of conventional 
orthotropic laminates, little is available which examines the postbuckling behavior of lam- 
inates like the [±50/35 ] 5 laminate used in the HiMAT effort. 
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1.2 Review of Previous Pertinent Work 

Considerable literature is available which documents analytical and experimental work 
involving various subsets of the general anisotropic plate with bending-extensional cou- 
pling. Leissa has presented a collection of such work [5], both analytical and experimental, 
and so it will not be discussed here. 

1.2.1 Global Approximation Function Techniques 

First, the work of Ashton and Whitney in the early 1970’s is of interest because these 
two authors both examined the analysis of anisotropic plates (in the bending sense) from 
the point of view of global approximation functions. Ashton took the classical Rayleigh- 
Ritz approach, while Whitney worked with a modified Fourier series expansion technique. 
In both cases, the determining factor in the accuracy of the results was whether or not the 
natural boundary conditions were satisfied. 

Ashton’s work [6,7] is limited to the analysis of a few selected problems whose ge- 
ometric boundary conditions can be satisfied using simple trigonometric functions. In 
the problems where natural boundary conditions are satisfied by the chosen approxima- 
tion functions, convergence to accurate solutions is good; however, in those cases where 
natural boundary conditions are not satisfied, convergence is shown to be unsatisfactory 
even for gross response quantities such as maximum displacement and buckling load [8]. 
This limitation is severe in problems where the choice of separable functions of the form 
X(x)Y(y) satisfying natural boundary conditions is precluded by coupling between bend- 
ing and twisting. If coupling between bending and extension is also present, satisfaction 
of natural boundary conditions becomes even more difficult. 

Later work by Whitney [9,10] shows that the use of a modified double Fourier series 
analysis results in better convergence because functions can be derived which satisfy natural 
boundary conditions. For the problems studied, the technique works well; however, the 
method is still not generally applicable in an automated fashion to a wide range of problems, 
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since new functions must be derived for each set of boundary conditions. Also, if bending- 
extonsional coupling is present, the method as it is presented is not applicable. 

The use of the exterior penalty function provides an easily automated way to use global 
approximation functions for the analysis of a wide range of problems involving complex 
mechanical couplings and boundary conditions. The penalty function method allows the 
choice of any set of approximation functions, so long as they are separable and adhere to 
simple symmetry rules. Boundary conditions are satisfied automatically, freeing the analyst 
to select approximation functions which represent complicated response quantities more 
accurately than would functions limited by boundary condition requirements. This freedom 
from the limitations of boundary conditions should make possible the use of approximation 
functions chosen for accuracy, making the Rayleigh-Ritz technique a viable alternative to 
the finite element method for problems having simple shapes. 

1.2.2 Reduced Basis Work 

Considerable work has been devoted to the development of different reduced basis anal- 
ysis techniques, although this work goes by many different names. Most of the techniques 
are methods by which large, computationally expensive equation systems describing a dis- 
cretized structure are approximated by much smaller sets of equations. The techniques are 
generally based on the approximation of the problem solution by a combination (whether 
linear or nonlinear) of a set of basis vectors assumed to span the desired solution space. 

The different methods are distinguished primarily by the differences in the basis vec- 
tors themselves. The modal techniques used in the structural dynamics field are a form 
of reduced basis analysis, since the problem unknowns are being approximated by a linear 
combination of modes which are typically selected to be the eigenvectors of the problem 
at hand [11]. Fairly recently, basis vectors have been developed for use in static nonlin- 
ear analysis; these vectors often consist of derivatives of response quantities with respect 
to some generalized path parameter [12]. More recently, Noor has analyzed anisotropic 
plates using various order derivatives of response quantities with respect to the anisotropic 
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material stiffness coefficients as basis vectors [18]. Elements of this last approach will be 
applied herein. 

Noor’s work in the area of anisotropic analysis has dealt exclusively with the finite 
element method, used in conjunction with the symmetry relations published in reference 
[13]. Taking this approach, he has simplified the process of obtaining a set of basis vectors 
for an anisotropic problem in three ways. First, the equation systems defining the basis 
vectors are shown to be approximately half the size of the full equation system; second, 
only one set of nonlinear equations must be solved to obtain the basis vectors. Finally, 
many elements of the reduced system are shown to be zero for symmetric or antisymmetric 
loading. 

However, the finite element method is a technique which is best suited to complex 
geometries not easily characterized analytically. For geometrically simpler structures, a 
classical Rayleigh-Ritz analysis can take the place of a finite element model. It is partic- 
ularly interesting to consider such an analysis in the context of the symmetry relations 
detailed in [13] since these relations can be directly applied in the choice of approximation 
functions. These simplifications find use in all stages of the analysis, from formation of the 
full and reduced equation systems, to appropriate choice of approximation functions. 
1 . 2.3 Experimental Work 

Jensen [14] studied the buckling and postbuckling behavior of fully anisotropic plates 
with bending-extensional coupling in an experimental context. Square plates of various 
laminations were tested, and the results compared with both classical Rayleigh-Ritz analy- 
sis and finite element analysis results. (Again, the Ritz functions met geometric boundary 
conditions.) The eleven layup configurations examined ranged from orthotropic to fully 
anisotropic with bending-extensional coupling. 

The work is fairly extensive, but the results do leave room for additional investigation, 
as noted by the author in his recommendations. The laminates studied in reference [14] 
were carefully chosen in order isolate several different distinguishing characteristics inherent 
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in panels with anisotropy and bending-extensional coupling. With this groundwork laid, 
it is now desirable to study laminates which might be used in the design of actual aircraft. 
It would also be desirable to test anisotropic laminates using a test setup which has been 
employed extensively by others (for example, [15]) for two reasons. First, the suitability of 
such test methods to the testing of anisotropic panels with bending-extensional coupling 
could be evaluated; second, test results could be qualitatively compared with the results 
of others. 

1.3 Objectives and Scope 

The overall objective of the present work is to study the effects of anisotropy and 
bending-extensional coupling on the buckling and postbuckling response of flat panels. 
The specific objectives of the present work are: 

1. To develop a modified Rayleigh-Ritz analysis technique and apply it to the linear 
and nonlinear analysis of anisotropic panels with bending-extensional coupling. 
Modifications include, for the linear case: 

a. Exploiting known symmetries of anisotropic panels in the selection of ap- 
proximation functions [13]. 

b. Developing and applying a reduced basis technique based on these same 
symmetries. 

c. Enforcing the geometric boundary conditions via an exterior penalty function 
approach, rather than by choosing approximation functions which satisfy the 
boundary conditions automatically. 

For the nonlinear case, only modifications (a) and (c) will be used. 

2. To gain insight into the postbuckling response and failure characteristics of pan- 
els of various aspect ratios which possess anisotropy and bending-extensional 
coupling by conducting postbuckling experiments involving such panels. Any 
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observed phenomena which are not present in the behavior of panels without 
anisotropy or bending-extensional coupling will be noted. 

The scope of the current work includes: 

1. Analysis of thin, flat laminated composite panels possessing anisotropy and bend- 
ing-extensional coupling using both the modified Rayleigh-Ritz technique and 
readily available finite element analysis tools. The effects of transverse shear de- 
formation are included in the Rayleigh-Ritz analysis, but not in the finite element 
analysis. Failure criteria will not be considered. 

2. Uniaxial loading in all analyses and experiments, with transverse loading being 
used as a check case in the Rayleigh-Ritz analyses. 

3. Use of simple linear and trigonometric approximation functions in the Rayleigh- 
Ritz analysis. 
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Chapter 2 

Problems Involving Bending-Extensional Coupling and Anisotropy 


2.1 Mechanical Couplings 

The concepts of anisotropy and bending-extensional coupling are most easily explained 
in terms of the constitutive relations for a laminated panel. The following constitutive 
equation relates the extensional stress resultants N 1 , N 2 , and N 12 , bending resultants Mi, 
M 2 , and M 12 , and transverse shear stress resultants Q 1 and Q 2 to the corresponding strain 
quantities defined below via the material stiffness matrices A, B , D, and C: 
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(Note that for purposes of calculating the material stiffness coefficients shown above, the 
reference surface is taken to be the middle surface of the panel.) For later reference, the 
constitutive equation can be written in a more compact form: 

*> j = 1, 2, 6 

(2.2) 

Cfc/J \ £s J k, l = 4, 5 

For convenience, the matrix C is partitioned as shown. The A coefficients axe known as 
the extensional or membrane stiffnesses of the panel. The D coefficients are referred to as 
the bending stiffnesses, and the B coefficients are called the bending-extensional coupling 
stiffnesses since they couple bending and extensional behavior. Finally, the C terms are 
the transverse shear stiffness coefficients. All these coefficients depend upon the layup 
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configuration and are derived from basic theory of laminated composites as may be found 
in any of the standard texts on the subject (See, for example, [8].) 

A modified von Karman type nonlinear plate theory is used with transverse shear 


strains and bending-extensional coupling included. The strain-displacement relations are: 

Extensional Strains Curvatures 


£1 = diui + %(diw) 2 
e\ = d 2 u 2 + \{d 2 w ) 2 
2e?2 = $i«2 + d 2 ui + diwd 2 w 


kj = d\4>i 

k 2 = d 2 <f> 2 
2«i2 = d\<f> 2 d 2 (j)\ 


(2.3) 


Transverse Shear Strains 


2ei3 = <f>i + d\ w 
2e 2 3 = <f> 2 d 2 w 

where the symbol d{ indicates differentiation with respect to a The conventions which 
are used for the plate geometry, force and moment resultants, loads, and displacements 
are shown in Figure 2.1. 

Three simplified forms of equation 2.1 can be obtained by examining three major 


types of laminates: 

1. Midplane-symmetric orthotropic laminates 

2. Unsymmetric orthotropic laminates 

3. Anisotropic (nonorthotropic) laminates 

The first category, that of midplane-symmetric orthotropic layups, possesses a consti- 


tutive matrix of the following form: 


Ni 

N 2 

n 12 



M\ 2 


U:J 


An A12 

A-12 A22 


Symmetric 


D n 
D\ 2 


D\ 2 

d 22 

£>66 


'44 


C 55 


fc 2 

2^12 

< *1 * 
« 2 
2«i2 

2£23 

. 2£i3 , 


(2.4) 
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This form occurs because a midplane-symmetric laminate, by definition, possesses mirror 
symmetry with respect to the midplane, causing the entire B matrix to vanish. The only 
remaining coupling is embodied in the A\ 2 and D 12 terms. This coupling is known as the 
fami liar Poisson effect, i.e., the occurrence of lateral strain when an axial load is applied, 
and vice-versa, as well as the presence of curvature in both directions. 

When an orthotropic laminate is not midplane symmetric, the B,j terms are nonzero 
and the following relation is obtained: 



This relationship, which exhibits coupling between bending and extension, will occur when 
the layer orientations or thicknesses on one side of the laminate midplane are not identical 
to orientations or thicknesses on the other side. The presence of the J3,y terms means 
that when the panel is subjected to inplane forces or displacements, curvatures are ob- 
served. Likewise, when bending-type loads or curvatures are applied, extensional strains 
are observed. 

In postbuckling analysis, the effects of bending-extensional coupling must be kept in 
mind for the following reason. Consider a geometrically symmetric panel with symmetric 
boundary conditions and bending-extensional coupling. When a symmetric axial load is 
applied, the prebuckling deformation is symmetric. However, the postbuckled configu- 
ration can be such that precise satisfaction of the symmetry conditions described in the 
next section does not occur. The response is either symmetric or unsymmetric; it is never 
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antisymmetric. This fact was noted during the finite element analysis of three of the six 
panels tested and is discussed further in Chapter 5. 

Bending-extensional coupling causes one additional difficulty. It means that more in- 
formation is required in order to evaluate moment and force resultants. In problems where 
there is no bending-extensional coupling, if only extensional strains are known, then it is 
a simple matter to obtain force resultants just from the constitutive relations; however, in 
the presence of bending-extensional coupling, curvatures must also be known to determine 
inplane force resultants. Similarly, when bending is experienced, moment resultants are 
easily found using just curvatures and the constitutive relation. When bending-extensional 
coupling is present, extensional strains are also needed. This additional information re- 
quirement can be especially troublesome from an experimental viewpoint, since totally 
different techniques are used to measure curvatures and strains. An alternative to mea- 
surement of curvatures in addition to strains is the calculation of the curvatures and middle 
surface strains using the assumption that the strain variation is linear in the thickness di- 
rection: 

£ 11 ] pii , | f*n 

£22 > = < £22 '+M *22 

£l2 X3=z V £12 J l *12 

where 2 may take on any value from —h/2 to h/2. 

Anisotropic laminates can be distinguished from orthotropic laminates by the presence 
of any of the following four sets of coupling coefficients: 

1) Inplane shear-extensional coupling (Ai 6 , A 2 6 terms) 

2) Bending-twisting coupling (D ig, D 2 6 terms) 

3) Extension-twisting coupling (B 16 , B 2 t terms) 

4) Transverse shear coupling (C 45 term) 

Henceforth, these material stiffness coefficients will be referred to as the anisotropic stiffness 
coefficients, and are not to be confused with the bending-extensional coupling coefficients 
for orthotropic laminates, discussed above. 
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Inplane shear-extensional coupling causes the panel to experience inplane shear defor- 
mations when it is subjected to tension or compression. A schematic is shown in Figure 
2.2 which exaggerates such deformation for the case of a panel subjected to uniform axial 
loading. The panels analyzed in the present work are subjected to uniform end shortening, 
so that the loaded edges are not allowed to shear, as shown in the figure; however, the 
effects of shear-extensional coupling are still observed. They are particularly obvious in 
the postbuckling range; the nodal lines of the buckle pattern tend to skew diagonally, as 
shown in the contour plot of Figure 2.3. 

Bending-twisting effects become obvious when the plate begins to exhibit out-of-plane 
behavior; the plate will then become twisted (nonzero K 12 ). During a test, this behavior is 
easily observed by placing displacement transducers at geometrically symmetric locations 
and observing the difference between the two displacements. (See Figure 2.4.) This type 
of coupling also contributes to the overall skewing behavior shown in Figure 2.2. 

Extension-twisting coupling can occur only in laminates which axe not midplane 
symmetric, since it is a form of bending-extensional coupling, and, as mentioned above, 
bending-extensional coupling occurs only when the layup is not midplane symmetric. The 
result of this type of coupling is that the plate experiences out-of-plane warping or twisting 
displacements when it is subjected to an inplane load. In the lab, this behavior can be 
observed in the same way as bending-twisting coupling; however, if both types of coupling 
are present, their effects cannot be isolated from each other in this way. 

The last type of coupling exhibited by anisotropic panels is transverse shear coupling. 
This coupling, unlike all the others, does not have an obvious, physically observable ef- 
fect on the panel response. However, it has recently been demonstrated (see [16]) that 
transverse shear can have a pronounced effect on the postbuckling response of laminated 
composite panels. Therefore, it is expected that coupling between the two transverse shear 
quantities can also be important in determining the response. 
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2.2 Symmetries Exhibited by Anisotropic Panels 

The symmetry of the response of an isotropic or orthotropic panel is simple: if the 
panel geometry, loading and boundary conditions have mirror symmetry, then the entire 
panel response, including displacements, stresses, strains, and force resultants, will exhibit 
the same mirror symmetry. Similarly, if the loading is perfectly antisymmetric, then the 
response is also antisymmetric. However, for anisotropic panels, the situation is not so 
simple; inversion symmetry rather than mirror symmetry is exhibited. 

Inversion symmetry, according to [13] and [17], means that if geometry, loading, and 
boundary conditions are symmetric, then the displacement field will satisfy the following 
symmetry statements: 

i,x 2 ) = -u a (-xi,-x 2 ) 

w(xi,x 2 ) = w(—x\, —x 2 ) a = l,2 (2.7) 

^o(®l>®2) = ^e*( *1 » ®2) 

Symmetry relations for the stress resultants are: 

N a p(x 1 ,X 2 ) = Naf)( 3?l, x 2 ) 

M a p(xi,x 2 ) = M a p(—xi,—x 2 ) a, P = 1,2 (2.8) 

Qa/}( x 1 > •*-2) — Qotfii, 3-1? * 2 ) 

This kind of symmetry is satisfied for both orthotropic and anisotropic panels. In a post- 
buckling problem where the out of plane response may be antisymmetric, the same sym- 
metry conditions are obeyed, except that 

w(xi,x 2 ) = —w(—x i,—x 2 ) 

lj®2) = ^or( ®1) ^ 2 ) 

This statement is only true for the case of no bending-extensional coupling. When bending- 
extensional coupling is present, there are no straightforward symmetry conditions unless 
the buckle pattern is symmetric (i.e., w is symmetric). 
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Before considering exploiting the aforementioned symmetries in the classical Rayleigh- 
Ritz technique, one must understand the difference between symmetries exhibited by or- 
thotropic and anisotropic panels. This is most easily done by considering how the problem 
would be solved when using the Rayleigh-Ritz method. Consider how each displacement 
degree of freedom might be approximated. Notice that the symmetry conditions on u a 
and 4> a axe totally satisfied if these displacements are approximated by even functions in 
one of the coordinate directions and odd functions in the other. In fact, in the solution of 
orthotropic problems, these are precisely the kinds of functions which are used. The sym- 
metry of the problem requires association of function symmetry with coordinate direction; 

\ 

for symmetric loading and therefore symmetric w, the symmetry conditions are: 


r Uy ' 


r -«i ' 


' Ul ' 

u 2 


U 2 


-«2 

W 

= i 

W 

• = < 

W 

<j> 1 


1 


4>i 

■ 4> 2 ‘ 

(*1,Z2) 

< 4> 2 - 

(-Xl,Z 2 ) 

■ —<t>2 > 


(2.9) 

where the subscripts on each vector denote evaluation of the vector at the point shown. 
If the problem, and therefore w, is antisymmetric and there is no bending-extensional 
coupling, then 


f «i ' 


’ Ui 'j 

■ -It! > 

I U2 


-«2 


■ W 

* = < 

—W > = < 

—IV 

4>\ 


h | 

-h 

l <f) 2 > 

(*1,*2) 

k -<^2 ) (- Xl , X2 ) 

<f>2 ) 




( 2 . 10 ) 


These symmetry relations are the same as the modal symmetry relations cited in [18]. 
Using these symmetry relations, convenient tables can be formed indicating how the ap- 
proximation functions should be chosen: 
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Table 2.1. Orthotropic Panel Symmetry States 



Symmetric w 

Antisymmetric w 

Xl 

X 2 

X\ 

X2 

w 

even 

even 

odd 

odd 

u\ and (j > i 

odd 

even 

even 

odd 

t /2 and <f ) 2 

even 

odd 

odd 

even 


Henceforth, the two distinct sets of functions will be referred to as corresponding to 
an “even state”, where w is symmetric, or an “odd state”, where w is antisymmetric. 

In the anisotropic case, the displacements still satisfy the inversion symmetry condi- 
tions, but they can no longer be approximated by a set of strictly even or strictly odd 
state functions as delineated for the orthotropic case. Rather, they must be represented 
by a combination of even and odd function sets since the two states become coupled. This 
coupling occurs due to the anisotropic material stiffness coefficients; when these terms are 
zero, the response is decoupled into a symmetric response due to a symmetric load, and an 
antisymmetric response due to an antisymmetric load. Therefore, as might be expected, 
the orthotropic problem is a special case of the anisotropic problem as far as symmetries 
are concerned. This fact will be used in Chapter 3 when operator splitting is introduced 
as a modification of the classical Rayleigh-Ritz method. 

Another consequence of the more complex symmetries exhibited by anisotropic panels 
is that the finite element model must be larger. For an orthotropic panel, only one quarter 
of the panel must be analyzed, since the two panel centerlines are lines of response sym- 
metry. However, for an anisotropic panel, the panel diagonal is the only line of response 
symmetry, and so half the panel must be analyzed. Reference [17] outlines procedures 
which can be used to analyze half of an anisotropic panel. Furthermore, as noted in 
later chapters, accurate representation of the more complex response of anisotropic panels 
requires a higher degree of discretization than for an orthotropic panel. 
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2.3 Governing Differential Equations of an Anisotropic Plate 
And Bending-Extensional Coupling 

A short discussion of the governing differential equations associated with the problems 
solved herein is presented here. The methods used throughout this thesis are energy meth- 
ods, but it is also useful to look briefly at the differential equations, since it is frequently 
helpful to view most problems from different perspectives in order to better understand 
and solve them. 

The static differentia] equation of an orthotropic plate with no bending-extensional 
coupling subjected to transverse load is solvable for a fairly large number of cases. Only 
even-order derivatives of displacements are present in this equation; it can therefore be 
solved analytically by the method of separation of variables. To effect such a solution, 
transverse displacement w is represented by the product of two functions, F(xi)G(x 2), 
where one of either F or G must be a trigonometric function, i.e., a function whose second 
derivative is a constant times the original function. 

However, when bending- twisting terms are included (i.e., terms involving Die and 
£>26), the equation becomes much more complicated due to the presence of odd as well 
as even-ordered derivatives. In this case, separation of variables will only be effective 
if either of F or G is of the form ke' mx , where k and m are constants and x is either 
x\ or x 2 . Even so, the resulting ordinary differential equation is difficult to solve, since it 
contains imaginary terms. The only known analytical solution to a plate problem involving 
bending-twisting coupling involves a clamped, elliptical panel under uniform load, which 
is certainly a very limited case. 

For an orthotropic plate subjected to inplane loading, there are three uncoupled differ- 
ential equations. Two of these equations are simply the plane stress equilibrium equations 
and are solved first to give the inplane displacement field. The third equation involves 
only the out-of-plane response and is identically zero if there is no transverse load. How- 
ever, if there is bending-extensional coupling, then the coupled constitutive relationship 
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serves to couple all three equations, and the solution becomes commensurately more diffi- 
cult. If anisotropy is included, the additional difficulty of the odd-order derivatives again 
presents itself. Of course, if the nonlinear terms axe added to the strain-displacement re- 
lations, then the inplane and out-of-plane differential equations are coupled even without 
bending-extensional coupling, since the force resultants will then contain terms involving 
w. 

Therefore, the differential equations of an anisotropic plate with bending-extensional 
coupling are sufficiently complex to warrant numerical solution. Available techniques in- 
clude Galerkin’s method and, from an energy standpoint, either the Rayleigh-Ritz tech- 
nique or its discretized version, the finite element method. Only the two energy methods 
will be considered herein. 
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Chapter 3 

Development of the Modified Rayleigh-Ritz Technique 


3.1 Review of the Rayleigh-Ritz Technique 

The Rayleigh-Ritz method is a well-known technique for approximating the behavior 
of engineering structures. It is the forerunner of the finite element method, in that it 
treats the entire structure as one large finite element. The classical method begins with 
the potential energy description of the problem, and results in a set of equations to be 
solved for the amplitudes of chosen approximation functions. This section reviews the 
technique and introduces the notation which is used in succeeding sections. 

3.1.1 Development of the Exact Energy Expression 

The first step in the solution of the plate problem which is the topic of this thesis is to 
write an expression for the energy of the plate, where W ext is the work of external forces: 

U = i [ {Oijeij + cr j3£*3 + a 3i e 3i )dV - W* xt i,j = 1, ... ,2 (3.1) 

* Jv 

Using the strain-displacement relations shown in equation 2.3, the plate strains can be 
written as the sum of a linear and a nonlinear part as follows: 


Linear terms: 


' £u 

L 

f £ 11 I 

L 

*11 


d\U\ + x 3 di<j>i 

£22 


^22 


K 22 


dzU2 + X 3 d 2 (f>2 

2£i2 

* = < 

o-O 

Z£ i2 

* + X 3 < 

2«12 

> = < 

d\U2 + + X 3 (di4>2 + d2<i>\) 

2^13 


2^13 


- 


<t>i 4- d\ w 

. 2^23 , 


. 2£23 , 


k > 


„ <^2 + @2™ J 


(3.2) 


Nonlinear terms: 


' £ 11 

NL 

f £?i 1 

NL 

r } 

622 

r 

i e ° 2 i 

. = i 

l&w) 2 

„ 2£i2 > 

1 

l 2£? 2 J 


{ d\wd 2 w J 


(3.3) 
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Substituting these expressions into the strain energy and integrating over £3 introduces 


force and moment resultants into the strain energy expression: 


V = \t 

( 

( diui 

[JV, N 2 1V 12 ] 

d 2 u 2 

2 J 

V 1 

[ diu 2 + d 2 ui t 


f ) 

+ [Mi M2 M12] < d 2 <j > 2 / 

[di<l >2 + d 2 (j>i J 


+ [Ni 


N 2 


^ 12 ] < 


i(diw) 2 

§(0 2 u,) 2 


[ d\wd 2 w 


(3-4) 


+ ^{VX}) dA - w 

The following vectors axe now defined for ease of notation: 


' dim ' 


\d\w ) 

d 2 u 2 


\d\w 

diu 2 + d 2 v.i 


di wd 2 w 

di<f> 1 
d 2 <f> 2 

{£ NL } = 1 


d 2 <f>i + di<f> 2 


■ 

<j>i + ditv 


• 

„ <i> 2 + d 2 w , 


K * * 


These definitions transform the strain energy into a more manageable form: 


(3.5) 


U 


(3.6) 


= \J {[ Ni 1 V 2 N12 I Mi M 2 M 12 I Qi Q 2 ]{e L ) 

+ [Ni n 2 isr 12 ]{£ NL » dA - w ext 

The force and moment resultants axe related to the strains by the constitutive relations, 


which are: 


N l > 


/ c-0 \ 

n 2 


e 2 

n 12 


Oc- 0 

^12 

— - 


A B • * 


— 

Mi 

* = 

B t D ■ 

i 

Kl 

M 2 


• C 


«2 

Mi 2 


L J 


2ki 2 

Qi 


2e 2 3 

. Q 1 , 


. 2£i3 , 


[C]{{ £ L } + {e NL }} (3.7) 


19 


This substitution finally leads to the following concise expression: 

V = \ J [[ £ l 1 + i£ NL ]] [C] { {e L } + {« NL }} dA - W‘" 

= | / [[e L ] ( C]{ £ l } + [ £ L ][C]{ e NL } (3.8) 

+ [e NL ][C]{e L } + [e NL ] [C] {e NL }] dA - W ext 

The energy is now composed of four distinct parts: terms involving either {e L } or 
{e NL } only, or combinations of {e L } and {e NL }, and the work of external forces. Now the 
energy can be conveniently subdivided in the following way: 

u - u 2 + u 3 + u 4 - w ext 

where 

U2 = lJ[e L UC){e L }dA 

V» = \J \ (^ L ] [Cl {e NL } + [e NL ] (C) {e 1 -} dA (3.9) 

U 4 = \ J [£ NL ] [C] (£ nl ) dA 

The numeric subscripts indicate the power to which the generalized displacements are 
raised in each of the three terms. Notice that, due to the definitions in equation 3.5, the 
strain energy is now written entirely in terms of the material stiffness coefficients and the 
five displacement degrees of freedom u\, u 2 , <f>i, <j> 2 , and w. 

3.1.2 Selection and Use of Approximation Functions 

At this point, the approximation functions for each of the five displacement degrees 
of freedom are inserted into the energy expression. In the classical version of the Rayleigh- 
Ritz method, these functions must be: 

1. Capable of satisfying all geometric boundary conditions 

2. Linearly independent 

3. A subset of a complete set of functions (e.g., Fourier series) 
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4. Differentiable up to the order present in the strain energy 

The chosen set of approximation functions can be written as shown in equation 3.10 
below, where X is a vector of as yet unknown coefficients of all the approximation functions, 
and F is a vector of the functions themselves. The indices P, Q, R, S, and T range over 
the number of functions chosen to approximate each corresponding displacement; all the 
ranges may or may not be the same. (After the array multiplication, all repeated subscripts 
denote summation.) 


' «i 1 
«2 


*1 \ =[F]{X}= [F? Fp Fp Fp 
<t>2 

K W J 


X“ l 

xp 

x? 


I x\ 


(3.10) 


Substitution of these functions into equation 3.5 yields approximate expressions for 
{e L } and {e NL }. The linear strains then become 


{£ L } = J&Xj = 


■d 1 Fp 

diFp 
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d*Fp 

diFp 
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d\Fp 
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• 

d 2 Fp 


• 


d 2 Fp 

SiFp 
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Fp 
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d x F$ 

_ 


• 

Fp 

d 2 F$ 


(Xp 

X$ 

Xp 

Xp 

l x% 


(3.11) 


where m is the sum of all the ranges of P , Q, R, S , and T, i.e., the total number of approx- 
imation functions used, and a ranges from one to eight, the number of strain quantities 
being used herein. The nonlinear strains are slightly more complicated. Approximation 
results in the three dimensional array rather than a simple two-dimensional matrix. 

This array can be thought of as having eight “planes”, where each plane corresponds to 
the index a taking on a value ranging between one and eight. Therefore, each “plane” 
of this array corresponds to one element in the nonlinear strain vector shown in equation 
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3.5. There are only three nonzero planes, since there are only three nonzero nonlinear 
strain terms. Each plane is a matrix consisting of products of appropriately differentiated 
approximation functions; for example, the first plane, corresponding to the strain term 
= \{d\ u>) 2 is as follows: 



■Xp ' 

T 

. . . . 


Xp' 1 


X P 




Xp 

[eft] = X/ftXj = . 

X*' 

' 


< 

xp 


xp 




xp 





rn,m 

Vll> 

a t 2 


(3.12) 

where T\ and T2 have the same range as T, i.e., the number of functions approximating w. 
All the zero rows and columns are included to maintain compatibility with equation 3.11, 
so that the indices I and J still have the same ranges. 

At this point, it becomes more convenient to abandon the matrix notation of equations 
3.11 and 3.12 in favor of index notation, so that reference is made to the two R arrays 
mentioned above, rather than to the detailed contents of F. Using this notation, the two 
strain vectors become 

{e L } = R L aI X T 


(3.13) 

{e NL } =^L jX/Xj 

The matrix of material stiffness coefficients is then denoted by C a 6, where b, like a, ranges 
from one to eight. (Recall that when index notation is used, repeated indices denote 
summation. In the present work, this rule applies only to subscripts.) Furthermore, the 
arrays K, F, and G may be introduced as the results of the integration over the area of 
the panel, excluding the unknowns X/. With these modifications, equation 3.9 becomes: 
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u = u 2 + u 3 + u 4 - w ext 

= fjX/Xj + -FukXiXjXk + -GijklXiXjXkXl - W ext 

where 


U 2 = \ JxjR^C^RljXjdA 
= ±K IJ X I Xj 

U* = \ [ {XjR^CatR^XjXK + XjXjR^CatR^X^dA 

_ (3-14) 

= 2 FukXjXjXk 

lh = \j XjXjR^ a C ab R^ L X K X L dA 

= \g ijkl X i XjX k X l 


3.1.3 Generation of the Nonlinear Equations in Xr And the Recursion 
Formulas for the Newton-Raphson Iterative Procedure 

The next step is to take the variation of the energy with respect to the unknown 

coefficients and require that the coefficient of the resulting variational quantity 6 Xj vanish 

in order to satisfy the requirement for minimum potential energy. That is, 


6U = —^-SX! — 0 




dU 

dXi 


= 0 


(3.15) 


A linear set of equations K jjX j — Pi is obtained if the higher order U$ and U 4 terms 
are omitted from the energy; otherwise, a nonlinear set of equations in X/ is obtained. A 
typical method of solving such a system of equations is the well-known Newton-Raphson 
technique, described briefly in Section 4.5. When using the Newton-Raphson technique, 
one solves a system of linear equations in each iteration. The equations have the form 


9 f dU \] AV / r> 9U 

ax,(axJj A J { P/ ax, 


(3.16) 


Thus, in order to solve the problem, two partial derivatives of the energy must be taken. 
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3.1.4 Svmmetrization of the Nonlinear Arrays 


To generate both the linear and nonlinear Newton-Raplison equations correctly re- 
quires some care with the index notation. Appendix A describes briefly the correct proce- 
dure, referred to as symmetrization, using only the quadratic contribution to the energy. 

If the unbarred symbols K , F, and G are now used to indicate symmetrized arrays 
as defined in Appendix A, the following nonlinear equations are obtained, along with the 
corresponding Newton-Raphson recursion formula in X/, where Pi is the load vector: 

Nonlinear equations: 

KuXj + ~FukXjX k + \g IJK lXjX k X l - Pj = 0 (3.17) 

Newton-Raphson recursion formula for iteration r + 1: 

[K,., + f,j K x’k + g,jk tXjrxqAxy” 

= {P, - K,jX r j - F IJK X r jX r K - G UKL X r jX' K Xl} 

(3.18) 

3.1.5 Generation of Load Vector P 

The load vector P was generated from the variation of the work of external forces 
W ext and is therefore dependent upon the particular load system being considered. For 
the case of a uniform transverse load, the work done is particularly simple; it is simply 
the product of the load magnitude and the transverse displacement, integrated over the 
loaded surface: 

pwdA = J P X™Ft dA (3.19) 

Then, the load vector P shown in equation 3.17 is just the variation of pF ext : 

SW exi = P T 6Xf =p J FfdA 8X% (3.20) 

In the case of an axial load, the situation is only slightly more complex. The work of 
an axial load is the product of the load magnitude and the axial displacement, integrated 
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along the loaded edge. If it is assumed that the panel is loaded along the edges x\ = ±|. 
then the work is written 


W ext = f Ni ui(-^,x 2 ) + Ni ui(^,x 2 )dx 2 

£ 

= r ' N t X'p' [f?(-|,*,)+ F?(|,* a )] d*: 

J 2 

Finally, taking the variation results in the load vector P: 


(3.21) 


6W 


ext 


= PpSXp = J^Ni [f^(-|,* 2 )+ F^(~,x 2 )] dx 2 6Xp 


(3.22) 


3.2 Enforcement of Boundary Conditions Using 
A Penalty Function Approach 

In section 3.1.2, it was pointed out that in the classical version of the Rayleigh-Ritz 
technique, all the approximation functions must satisfy the geometric boundary conditions 
of the problem at hand. This requirement can be quite limiting, in that first, it is sometimes 
very difficult to devise functions which satisfy all the geometric boundary conditions, and 
second, those functions may or may not possess good convergence properties. Therefore, 
at this point the first deviation from the classical theory is made in the following two ways: 

1. Choose functions which do not necessarily satisfy any boundary conditions 

2. Enforce the boundary conditions using a penalty function 

The penalty function approach consists of adding a penalty term P to the potential 
energy expression and then working with the variation of the augmented energy in deriving 
equations in X/; that is, 

U p — U + U 
SU P = SU + 611 = 0 ' 

The penalty term represents the degree to which a given constraint (in this case, a boundary 
condition) is not satisfied. The constraint can be written as 

n = rh(X 7 ) = 0 (3.24) 
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where r is a penalty parameter and h(Xj) is a function of the actual constraint. It is easy 
to see that if the penalty function is formulated such that it is always positive, then non- 
satisfaction of the constraint will increase the value of the augmented energy. Therefore, 
since setting the first variation of U p to zero and solving for X / constitutes minimization 
of Up , the penalty term is forced to zero and the constraint is satisfied. Note that the 
addition of the external penalty function can be viewed as the placement of a very stiff 
spring at points where constraints are to be imposed. The penalty parameter r is then 
analogous to the spring stiffness. 

Careful selection of the penalty parameter r is very important, as noted in reference 
[19] and confirmed by results presented later. Proper choice of r strikes a balance between 
two requirements: the penalty parameter must be large enough to force satisfaction of 
the constraints, and small enough not to numerically dominate the problem and cause 
significant loss of precision. 

As an example of constraint formulation, take the case of a w — 0 boundary condition 
along an entire edge Xi = 0. An integral constraint for this boundary condition may be 
written as follows: 


= / [w(0,x 2 )f 
Jo 


dxi 


(3.25) 


The integrand is raised to a power to ensure that the penalty function has a slope of zero 
when the constraint is satisfied, and so no discontinuity is introduced at the constraint 
boundary. Obviously, any power of w will satisfy this requirement; however, the use of 
an even power ensures that the constraint is positive if it is unsatisfied. Furthermore, an 

advantage of using the square of the constraint is that its variation is easily added to the 

/ 

linear stiffness matrix K, as will be shown subsequently. 

Noting that w is represented by a series of approximation functions as described in 
Section 3.1.2, the constraint may be rewritten as 


h = J F%(0,x 2 )F%(0,x 2 )X Tl X T2 dx 2 


(3.26) 
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where T\ and T 2 have the same range as T, i.e., the number of approximation functions 
representing w. If the penalty function is included in the derivation of the equilibrium 
equations as described above, the following term is then obtained: 


Sh = ^-«X Tl = j Fj?(0,z 2 )F£(0,* 2 )X Ta dx 2 8X Tl (3.27) 

This constraint can be recast in the following form so that it can be added directly to the 
linear term of the equilibrium equations: 


Sh = 



i 

x';: 1 
* 5 ’ 

X*' 

xp 



vw 

a t 2 

• • ■ [ F t, f t,] t ,t. 

"m,rrT 



(3.28) 


This integral-type constraint becomes particularly simple if the approximation functions 
are orthogonal; then, the matrix shown in equation 3.28 is diagonal, since the integral 
is zero when T) ^ T 2 . Notice that all that has been said about the edge constraint can 
be applied to constraints along a line (or even a curve) anywhere on the panel. Point 
constraints are even simpler, since they are not integrated; however, the constraint must 
still involve an even power since it must always be positive. 


3.3 Operator Splitting 

The operator splitting technique used herein allows two modifications to the classical 
Rayleigh-Ritz technique: 

1. Partitioning of the approximation functions into strictly even and strictly odd sets of 
functions. 

2. Splitting of the constitutive matrix into orthotropic and anisotropic parts. 

The details and consequences of these modifications will be explained in the current section. 
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3.3.1 Basic Concepts 


The justification for partitioning of the approximation functions into even and odd 
sets lies in the panel symmetry conditions outlined in Section 2.2. There it was noted that 
orthotropic panels exhibit purely symmetric or antisymmetric responses, depending on the 
symmetry of the loading and boundary conditions, while for anisotropic panels, coupling 
between the two types of responses occurs through the anisotropic (nonorthotropic) ma- 
terial stiffness coefficients. Taking these facts into account, it seems reasonable to choose 
functions for the anisotropic Rayleigh- Ritz analysis which are either strictly even or strictly 
odd, since it is known that the response can be represented by some combination of the 
two types of functions. Therefore, a partitioned set of approximation functions is defined, 
where each partition is of the same form as equation 3.10: 

F/X, = F?X/ + X" = [[r*]j I (3-29) 

The superscripts e and o refer to even or odd symmetry states, respectively. Actually, 
since indices i through / will correspond only to even functions, and indices m through 
q will refer to odd functions, the superscripts are redundant; however, they will enhance 
readability later on. The uppercase Latin indices range over both sets of functions, while 
the uppercase script letters 1 and M indicate the total number of even and odd state 
functions, respectively. Recall from Chapter 2 that the term “symmetry state” refers to a 
set of approximation functions chosen according to the following table, which is repeated 
here for convenience: 
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Table 2.1. Orthotropic Panel Symmetry States 



Symmetric tv 

Antisymmetric tv 

x x 

*2 

Xi 

x 2 

w 

even 

even 

odd 

odd 

u\ and <f>i 

odd 

even 

even 

odd 

U 2 and <f >2 

even 

odd 

odd 

even 


At this point, a new convention is introduced for the partitioned arrays. It should 
be understood that whenever two arrays having indices with different ranges are added 
together, as in 

X I = Xf+X° m (3.30) 


each array is to be expanded so as to be conformable for addition. For example, this 
equation could be rewritten using matrix notation as 


{X} w = 



) i 

f {0}* 

1 + 

— 

Jw 1 

, {X°}m 4 


(3.31) 


Utilization of this convention will result in more concise equations, since the zero arrays 
will not be explicitly written. 

It was also stated in Section 2.2 that coupling of even and odd responses occurs 
due to the presence of nonorthotropic material stiffness terms. In order to more easily 
understand and use this coupling effect, the problem is separated into two parts: an 
orthotropic part and a nonorthotropic part. This separation is easily accomplished by 
writing the constitutive matrix as follows: 


C a4 = C°t + 


(3.32) 


where the superscripts 0 and A correspond to orthotropic and anisotropic (nonorthotropic) 
material stiffness coefficients, respectively. (Recall that the term anisotropic, as defined in 
Chapter 2, refers to all the constitutive terms with subscripts “16” and “26”, as well as 
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the C 45 term.) In this equation, A is included as a tracing parameter for the anisotropic 
terms; that is, when those terms are to be included, A is set to one. Otherwise, A is zero. 
As a result of splitting the constitutive matrix, the procedures outlined in section 3.1 will 
lead to a set of nonlinear equations of the following form: 

K?jXj + i F° jk XjX,< + \tf JKL XjX K X L 

+ X(K^jXj+P^ jk XjX k + ^ jkl XjX k X l ) = P, 

(3.33) 

The combination of the special choice of approximation functions and the splitting 
of the constitutive matrix results in changes to the equilibrium equations which help to 
reveal the effects of anisotropy. To illustrate how these changes come about, the linear 
equations will be examined. Consider first how the linear stiffness array K ij is formed. 
The orthotropic part of the K array consists of the sum of the following four terms: 

f ( nL(e) t ~iO pL(e) . pL(e) r ~\0 pL(o) 

K IJ = I \ K ia ^ ab K bj + K ia ^ab^bn 

1 pL(o) /~iO pL(e) 1 pL(o) pO pL(o)\ 

+ U ma ^ab K bj + K ma ^ ab K bn ) 

where the parenthetical superscripts indicate the types of functions used to form the R 
arrays. The anisotropic part of the K array is similar, except that C ° 6 is replaced by C^ b . 

When the integration shown above is performed, an interesting thing happens. For 
the orthotropic array, all those terms formed from a mixture of even and odd functions 
vanish. In the anisotropic array just the opposite occurs; all terms which do not contain a 
mixture of terms drop out. Therefore, different combinations of even and odd superscripts 
also serve to identify whether the arrays are orthotropic or anisotropic. The two linear 
terms may then be written in a form analogous to equation 3.14: 


(3.34) 


dA 


K 


ij 


Ku 


= I («“'* c°, 


pL(e) 1 piAo; p 
R-hi A R mn ^ah-R 


L(o) (~\0 pL(oF 


b 1L bj T “ma ^'ab'Hn 


dA 




= / 


(3.35) 


C A pL (0) . nHOJ piA p. 

nbRh r, T- R ma l-' ab n b j 


‘ah lv, bn 


L(o) , 


?L( 


! e> ) 


dA 


= K e ° + K° e 

in 1 mj 
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The effect of the splitting process on the linear problem is more easily seen from the 


matrix form of the full linear set of equations, with the tracing constant included: 


K 


ee 

O' 


k: 


+ A 


X°n 


K 


mj 


ra -(gi »»> 


(After the array multiplication is carried out, repeated indices denote summation.) Here 
it is easy to see how the even and odd parts of the equations decouple for an orthotropic 
problem (i.e., where A = 0), giving a symmetric response for a symmetric load system and 
an antisymmetric response for an antisymmetric load system. Conversely, one can also 
see how the coupling between the symmetric and antisymmetric systems is accomplished 
when the anisotropic stiffness terms are present. 

From the above discussion, it is easy to make the incorrect assumption that whenever 
a mixture of even and odd functions is used to form the K, F , or G arrays, the result 
will be associated with the anisotropic elements of C a j,. Rather, one way to determine 
whether terms are orthotropic or anisotropic is to assign a negative one to each odd set of 
functions comprising a term and a positive one to each even set. Then, take the product 
of the resulting combination of positive and negative ones. If the result is positive, then 
the term is orthotropic; otherwise, it is anisotropic. For example, take the two F terms 
shown here: 

(1) 7%j = J C*, cj, < e) ) iA 

(2) F°‘l = j (*»;> C° C° #<">) iA 

Term (1) is anisotropic, since the combination ( eoe ) gives (1 • — 1 • 1) = —1, and term (2) 
is orthotropic, since ( oeo ) gives (—1 • 1 • —1) = 1. These separation properties arise from 
the form of the partial derivatives involved in the strain-displacement relations; if all the 
integrals are formed, it is easily seen that terms vanish in the various cases because the 
integrand is an odd function. This can be mathematically proven using group-theoretic 
methods; however, the proof is not presented here. Interested readers are referred to [20] 
for additional material on the subject. 
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3.3.2 Array Svmmetrization When Partitioning is Used 


As one might imagine, introduction of the partitioning scheme detailed in the previous 
section complicates derivation of the equations of equilibrium and the Newton-Raphson 
equations. Some additional bookkeeping is required in forming the various arrays, since 
there are now four K arrays, eight F arrays, and sixteen G arrays. (These numbers can be 
reduced by noting that some of the arrays are identical.) In addition, the symmetrization 
described in Appendix A becomes somewhat more complicated. To illustrate the differ- 
ences between partitioned and non-partitioned symmetrization, the linear portion of the 
problem will be re-examined in this section. The following presentation will be analogous 
to that of Appendix A, so that comparisons are easily made. 

First, the energy term U 2 , which involves only AT, is modified. This term must now 
be written to include the four partitions of AT, as follows: 


U 2 = - 


K* e jXfx ; + k\ 


00 X° 

mn^m^ n 


+ a(a: 


yc y~° 


+ KZjXtX 


V 


(3.37) 


Again, in order to take the variation of U 2 , dummy indices must be employed, so that 

SU 2 = 1^- 8X a + ^-6X fi (3.38) 


dX a a dX fi 

The index a has the same range as *, while has the same range as m. Note that since U 2 
has been partitioned, its variation is now the sum of two separate variational quantities. 
Now, the product rule of differentiation is applied, giving: 

1 


6U 2 =: 


— ee (OX 

Max 


« dXf\ (— 


— eo 8X\_ —oe 9Xj 

e n ' ■** mj^ rr, 


tn dx% 


m J m oye 

uyv a/ J 


sx* 


+ 


dX 


dX° 


^eo„.dXZ 


dXi 


KZ n + A + K mj ^xj 


6X°fi 


(3.39) 


After conversion of the partial derivatives to Kronecker deltas, this expression becomes 

eu, =i [ pni (Sic.x; + x's ja ) + a (k‘° 6, a x° n + x^xy*,)] sxi 


+ [JC. + xy„,) + a (JtZx:6„ + x“ sx$ 


(3.40) 
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Elimination of the Kronecker deltas gives 




1 r 


K“X' + K“x; + A (/C X° + KZ.XZ .) ] 
+ KZ f x ° m + a (t?,;x: + k';,x‘)] 


6X% 


sx; 


(3.41) 


(3.42) 


Finally, eliminating the dummy indices gives 

m =\ [((*“ + TT') x ; + A (lf“ + iC) -C) “7 

+ ((lC„ + iC.) xi + a (ir“ + 3C) *?) 

Regrouping terms and adding in the loading terms leads to the complete variational ex 
pression of the linear problem: 


SXf 


m =\ [(Ki + k'')x; + a (kZ + 7C,)x° m - p; 

+ \ [(3C„ + k°L)x°„ + mkZ + KZi)x; - k I ex. 


(3.43) 


Finally, because coefficients of the two variational quantities must vanish, the following 
two coupled sets of equations are obtained: 


\(K?j + K*J)x; + i A (KZ + KZi)X° m = p; 
\{K°Zn + KZ)X° + | a ( j ?“ + TZi)x; = P, 


(3.44) 


o 

m 


(3.45) 


With this step, symmetrization of the stiffness arrays K has been been achieved for a 
partitioned system. If the unbarred symbol K is now used to represent the symmetrized 
arrays, a set of partitioned equations identical to equation 3.36 is obtained: 

l yet \-e _i_ \ JZ'eo Y° D e 

+ AA j m A m x j 

Tf oo yo 1 \ jfoe ye no 

m 

Similar operations can be carried out on the F and G arrays in order to arrive at 
symmetrized arrays. For the same reason that there are only three unique K arrays 
(K^ = K^ t ), symmetrization results in only four unique F arrays and five unique G 
arrays. 
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3.4 Application of Reduced Basis Techniques in Conjunction 

With Operator Splitting 

Using the information presented in previous sections, one can obtain solutions to 
many different problems using classical linear and nonlinear equation solution techniques. 
However, additional information, as well as added computational efficiency, can be obtained 
by applying reduced basis techniques in conjunction with operator splitting. 

Implicit in a reduced basis technique is the assumption that the problem unknowns 
(in this case, displacements) can be approximated by some linear combination of a set of 
independent “basis vectors”. This assumption amounts to the transformation 

TO* = [rhMrWv (3.46) 

where each column of the matrix [F] is a preselected basis vector or assumed mode, the 
vector {A'} contains an amplitude for each mode, m is the total number of degrees of 
freedom contained in X, and Af is the number of basis vectors to be used. A common 
choice for the basis vectors has been derivatives of a nonlinear solution with respect to a 
generalized path parameter [21]; that is, 

= {if} {0} {#£}h.V (3-47) 

where r) is the path parameter and X is the nonlinear solution. These basis vectors axe 
commonly referred to as path derivatives. 

The transformation shown above is used to effect a reduction in problem size; thus, a 
much smaller system of nonlinear equations is obtained. For a finite element problem of 
some size, solution of reduced equations is of course much less expensive than solution of 
the full set of equations. The reduced equations are used in place of the full system until a 
pre-defined error measure dictates that the basis vectors must be regenerated using a new 
nonlinear solution. 
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Recently, work has been done in which the basis vectors are derivatives of a finite 
element solution with respect to a tracing constant A similar to that introduced in equation 
3.32 [18]. The same concept will be used here, so that the transformation matrix [T] is 
defined as 


PW 



m,Af 


(3.48) 


Note that each vector is evaluated at A = 0. This is because approximation of the problem 
unknowns in this way amounts to forming a “generalized” Taylor series about A = 0, 
in which the normally fixed coefficients of the various series terms are replaced by free 
parameters [21]. That is, the solution to the anisotropic problem is viewed as a large 
perturbation from the orthotropic solution. 

3.4.1 Generation of Basis Vectors 
3.4.1. a Linear Problem 

A set of recursive equations which define the basis vectors is obtained by successive 
differentiation of the equilibrium equations and evaluation of the result at A = 0, regardless 
of whether a linear or nonlinear problem is being solved. The linear problem will be 
discussed in the current section in order to introduce the concepts as simply as possible. 

First, it is more convenient to deal with one indicial equation, and so the two equations 
3.45 are combined to give a single equilibrium equation. 


KfjXf + I< mn 


Xi + A {K‘°Xl + K%X j) = P! + P° m 


(3.49) 


An equation defining the first basis vector, (X/) 0 = (Xf + A’^) 0 , is obtained by simply 
evaluating equation 3.49 at A = 0, giving 


(The zero subscript is a reminder that all the basis vectors are evaluated at A = 0.) Notice 
that this equation can be written as two sets of uncoupled equations in X * and X°. 

K'f (XJ) 0 = Pf 

TSOO t V°\ V>° 

A mnl A n/0 — m 

An equation defining the second basis vector is obtained by differentiating equation 
3.49 with respect to A one time and evaluating the result at A = 0. First, the differentiation: 


(3.51) 


„ dX\ 

Kff — ^ + K 


OO 

mn 


13 fiX mn dX 

Now, evaluation at A = 0 gives 


f)X° / fiX° dX* 

- + (KX + K™jXj) + A (K?: ^ 


L h 


0 (3.52) 


/ dX e \ / fix°\ 

\~ax) a + K ° m ° n \~d\)r (*;)„) (3.53) 


•> v sx 

Once again, an uncoupled set of equations is obtained: 

( d X]\ 

o = ~^°™ 0 


T, r oo ( dX \ T/’oe ( ve\ 

^ J o — ^mj\Xj) 0 


(3.54) 


These equations can be solved for the second basis vector, which is the first derivative of 
the problem unknowns with respect to A. Note that the right hand sides of the above 
equations contain (X°) Q and (JAj) o , which are known from the solution of equations 3.50. 

As many basis vectors as are needed can be obtained by continuing this process. The 
resulting equations are defined by the following recursion formulae: 


fd c X?\ 

o 




(d c X 
V dX 


■)- 


- c ( KJZ ( +K oe 


d c -'X e j 


dX c ~ 1 


fiX 


c = 1,2,..., jV 

(3.55) 
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Notice that the even and odd partitions will continue to be uncoupled in every set of 
equations. A large savings is factorization time is therefore obtained, since each partition is 
generally half the size of the full system of equations, and factorization time is proportional 
to the cube of the number of equations. Total factorization time is therefore reduced by a 
factor of four. 

An important simplification to the recursive equations can be realized if both the 
loading and boundary conditions are either purely symmetric or purely antisymmetric. 
Recall from Chapter 2 that the response symmetry of an orthotropic panel is the same as 
the symmetry of the load system; that is, if the load is symmetric, then the response is also 
symmetric, and similarly for an antisymmetric load case. Since equations 3.50 are simply 
the equilibrium equations for an orthotropic panel, pure symmetry or antisymmetry of the 
load system gives: 

symmetric loading =» P r " = 0 => (X£J 0 = 0 
antisymmetric loading => Pf = 0 =» (X®) 0 = 0 


These simplifications propagate through the recursive equations so that in each set of 
equations, only one matrix, either Kff or must be factored. For example, if the 

loading is symmetric, the recursive equations become 


Kli (*/)„ 

= Pt 

K 00 1 

( d c X°\ 

= -cK™i 


K d* Jo 

mj 

Kff 1 

(FXi\ 

= -cK£ 


K d\‘ ) 0 

in 


v-'Xf 


d 


C— 1 VO 


d\< 


/ 0 


c = 1,3,5, .. . 
c = 2,4, 6, . . . 


(3.56) 


The transformation matrix [T] also takes on a simpler form when the loading is sym- 
metric or antisymmetric. In the case of load symmetry, [T] becomes 


[r]m,A/- = 


{x e h ) 


o 




r Mil 
1 ax J 


( 3 . 57 ) 


M / 0 


m,Af 
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If the load is antisymmetric, then similarly, 


0 


[r]m,.V — 


((Whl 


1 0 } 0 ({w 1 } 


M ' 0 


(3.58) 


m,N 


3.4.1.b Nonlinear Problem 

For the nonlinear problem, the procedure for obtaining the recursion equations defin- 
ing the basis vectors is similar to that used for the linear problem. As with the linear 
problem, these recursive equations are only half as large as the original set of equations; 
furthermore, only the first of the equations is nonlinear. As before, equations defining 
the basis vectors are obtained via successive differentiation of the nonlinear equilibrium 
equations and evaluation of the result at A = 0. 

Recall first that the nonlinear equations obtained using the method outlined in Ap- 
pendix A are considerably longer than the linear equations, due to terms involving the 
three and four dimensional arrays F and G. In this section, only the terms involving K 
and F will be included in the explanation of the development, since addition of the G terms 
results in extremely long and unwieldy equations which do not significantly add to the un- 
derstanding of the concepts involved. The reader will find that patterns uncovered in the 
shorter equations presented here are very similar to those found in the longer equations. 

The procedure begins with the following coupled equilibrium equations, which already 
include the tracing constant A: 

rree v" c i / rpeee , rricee i T 7 ieee\ \re \re > jpooe i rioeo . rpeoo \ \ro \ro 
K ij ^j+K^ijk +- t jik + *kji) A >j A >k + Z \- t mni + *min + *imn) X rn A n 

+ A [KZX° m + (Fffc + F','° + FZ) + fJSi + 

= P‘ 


TF’OO -y 'O . ( -pooz » jpooe . jpoeo «_ rioeo i jpeoo , jpeoo \ ve y-o 

^mn^n ' rmni ' ” nmi ' ^ min ' r nim ' r imn ' r A n 

. \ \(T/‘oe. \re . / 77000 1 17000 1 77000 \yo yo i 0 ( l?eeo i t peoe ■ rpoee \ \re \re 1 

+ A [(A j + yr mn p + r n mp ' ** pnm/^n^p ' ijm ' * tmj ' ^ ^ j \ 


= PL 

m 


(3.59) 
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An equation defining the first basis vector, (X/) 0 = (Xf + is obtained by 

simply evaluating equations 3.59 at A = 0, giving 

Ktf (jq) 0 + {Ft’l + F'Z + Fff) (Xj) 0 (*{)„ 

+ 2 + iC?« + «) (^m)o TO„ = Pf 


r/'oo ( vo\ I ( TTiooe . jrtooe . 77 »oeo 

^ mn v-^n/O ' v-^mn* ' ^ nmi ' * min 


+ + *vr„ + fxj) W)« w)„ = a 


O 

m 


(3.60) 


(The zero subscript is a reminder that all the basis vectors are evaluated at A = 0.) These 
two equations are slightly more difficult to uncouple than the corresponding linear equa- 
tions were, in that assumptions must first be made about the load and boundary condition 
sy mm etry. Take the case of perfectly symmetric loading and boundary conditions, in 
which case the vector P £ is identically zero. To see the result of this, rewrite the second 
of equations 3.60 as follows: 


{K° m ° n + (F, 


ooe 

mm 


, jpooe , poeo 
i « nmt 1 min 


+ F, 


oeo , rrieoo , rpeoo / ve\ T>° 

nim + -Tjmn + * inm JoJ V A nio — r m 


(3.61) 


If P m is zero, then either the term in brackets is zero or (X°) 0 is zero. If (X£) 0 is zero, then 
the first of the recursive equations is just the equilibrium equation for the corresponding 
orthotropic problem: 

Kf,‘ W), + (Ftjt + P/a + FS/n W)„ (Xl)o = Pf (3.62) 

This equation completely defines the first basis vector. Note that it is nonlinear; however, 
it is also half the size of the original full system of equations. 
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A set of equations defining the second basis vector is obtained by differentiating 3.59 
with respect to A one time and evaluating the result at A = 0. First, the differentiation: 

aye / QX ? flye 

rsee i t TTteee i irieee i rrttc\ I j ve , ve k 

K ij ~gy + V*ijk + + *k> > > I ~g\ A >‘ +A > 


+ 2(F^ { + FXfn + F.Z°n) 
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/ TV'Oe V“ c 1 ( I pooo I T7I0OO I 77*000 \ Y”0 V° 

' l^mnp ' ^ nmp ' -^pnm j-^n^p 

+ 2 (^“ + FZ) + F° m %)Xt Xj) 


(3.63) 


Now, evaluation at A = 0 gives: 
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(3.64) 


ipoe (ye\ ( jpooo , rpooo , i pooo \ryo\ ( yo\ 

— ~&mi /0 — v x mnp + * nmp + * pnm) V A p/o 

- 2 (Fg“ + FZ) + FZ )) (X’X (X‘) 0 
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For a symmetric problem where (X °) 0 is zero, these equations become: 



/ f)X 0 \ 

y + + ^nmt + ^min + *nsm + * imn + )o ^ ~q^~ J ^ 

= - w)„ + 2(*g“ + jc; + w). (*;)„' 

(3.65) 

The first of these equations says that ^ ~^~J vanishes, while the second equation defines 

d ™ j . Notice that both equations are linear, since both (ATf) 0 
and (-X’m)o are known quantities. Also note that the same pattern of alternation between 
derivatives of even and odd parts of Xj is occurring, as was shown for the linear problem 
in equation 3.57. As before, as many basis vectors as are needed can be obtained by 
continuing this process. 

3.4.2 Generation and Solution of the Reduced Equation System 
3. 4. 2. a Linear Problem 

Once the desired basis vectors have been obtained, the reduction of the original system 
of equations can be done. Before beginning, it is convenient to rewrite the transformation 
matrix using the more compact index notation, as well as the superscripts e and o to 
indicate even and odd partitions as before. 

r [r'W 1 

[IW= =r;. + r^ e (3.66) 

- [^ 0 ]^,V J 

Using this notation, equation 3.46 can be written in greater detail as 

x! + x° m = (r* c + r° mc )x c (3.67) 

Application of the transformation is a fairly simple two step procedure. First, the 
displacement degrees of freedom in the original equilibrium equation 3.49 are replaced by 
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the transformation equation; then, the resulting equation is premultiplied by the transfor- 
mation matrix so that the resulting system is symmetric. 
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+ 1 di n in l nc 


+ rs 


dm 
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K 


(3.68) 


This equation is a symmetric system of Af equations whose solution gives the participation 
coefficients A c , and can be rewritten as 


K dc X c = V d (3.69) 

Once the participation coefficients are obtained, the full solution X j can be recovered using 
equation 3.46, and solution of the linear problem is complete. 

3.4. 2. b Nonlinear Problem 

The basis vectors generated by the procedure outlined in Section 3.4. l.b can be used 
to greatly reduce the computational effort involved in the static nonlinear analysis of 
an anisotropic panel. More precisely, the effort required for solution of an anisotropic 
problem can be reduced to slightly more than the effort required for the corresponding 
orthotropic problem, provided that the loading and boundary conditions are symmetric. 
The sequence of steps for the simplest possible application of the reduction technique to 
nonlinear problems is as follows: 

1. Given a load factor p, solve for the Af nonlinear basis vectors corresponding to 
that load level. Obtaining the basis vectors involves solution of one nonlinear 
set of equations (equation 3.62) and Af — 1 linear equations (similar to equation 
3.65). All these equation systems are half the size of the corresponding equations 
for the anisotropic structure. 

2. Compute the reduced nonlinear equations using these basis vectors. The resulting 
system will be much smaller than the original set of equations. Usually no more 
than fifteen basis vectors are needed. Also, since there is a definite pattern 
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defining the basis vectors which are zero, the computational effort required to 
generate the reduced equation system can be reduced by judicious use of this 
pattern. 

3. Solve the reduced system of nonlinear equations to get the participation coeffi- 
cients X c for each of the basis vectors. 

4. Recover the full anisotropic solution by taking Xj = Tj c X c . 

5. Choose the next load factor and repeat the process until the maximum desired 
load factor is attained. 

The nonlinear equation solutions delineated above can of course be performed using any 
suitable nonlinear equation solution technique, of which the classical Newton-Raphson 
method is an example. 

Computationally faster techniques have been developed which can be used in place of 
the above procedure; these techniques are outlined in [18]. There are two basic methods, 
one of which uses two successive single-parameter reductions, while the other applies only 
one two-parameter reduction. In both cases, one of the parameters is the anisotropic trac- 
ing constant used herein, and the other parameter is the load factor. Both techniques offer 
advantages over the modified classical Newton-Raphson technique described above, with 
the two-parameter reduction technique being the most efficient. The procedure described 
above has the advantage of requiring only one single-parameter reduction and is therefore 
somewhat easier to develop; from a computer programming standpoint, it should be used 
as a starting point from which to implement the techniques described in the references. 
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Chapter 4 

Numerical Studies of Modified Rayleigh-Ritz Technique 

4.1 Comments on Problem Selection 

The purpose of the present chapter is to present numerical results, both linear and 
nonlinear, obtained using the modified Rayleigh-Ritz technique of the previous chapter. 
Several different physical problems have been chosen for this purpose. In those problems 
involving anisotropic panels, the laminates, loading, and boundary conditions have been 
chosen to correspond as closely as possible to the two seven inch-wide experimental pan- 
els described in Chapter 5. Further comments regarding the detailed characteristics of 
each laminate may be found in Chapter 5. The only difference between the anisotropic 
panels presented in the present chapter and those of Chapter 5 is that in implementing 
the Ritz technique, some compromises had to be made in modeling the boundary con- 
ditions; these are discussed in greater detail later in the current chapter. As a result, 
quantitative comparison is not made between the Rayleigh-Ritz and experimental results; 
however, qualitative discussion is presented. Problems presented in this chapter which do 
not correspond to the panels of Chapter 5 are for comparison purposes only. 

4.2 Linear Results 

In Rayleigh-Ritz analysis, quality of results depends greatly on the suitability of the 
displacement approximation functions X,- to the problem at hand. The present section will 
determine the suitability of a very simple set of functions by presenting solutions to four 
problems involving various support /loading conditions and various layup configurations. 

The software written to implement the modified Rayleigh-Ritz technique allows the 
use of any set of functions satisfying the following two criteria: 

1. Separability: all functions must be of the form f(x)g(y). 

2. Symmetry: functions must be such that sets of functions can be formed which 
adhere to the symmetries shown in Table 2.1. 


44 



ORIGINAL PAGE IS 
K2B EGOR QUALITY 


The fairly simple set of functions shown in Table 4.1 satisfies these two requirements and 
is composed only of trigonometric functions such as sine and cosine, as well as simple 
linear polynomials. Various subsets of those functions were used to generate all the results 
presented in this section. 

To determine the suitability of these functions, two sets of symmetric support /loading 
conditions will be examined: 

1. Uniform transverse loading, with simple support boundary conditions. 

2. Uniform axial loading, with clamped/simple support boundary conditions. 
Schematics of these two problems may be found in Figure 4.1. 

Before proceeding, the term “simple support” for the case of a panel with bending- 
extensional coupling needs clarification. The classical definition of a simply supported 
boundary is that transverse displacement ( w ) and normal moment (M n ) must vanish along 
the boundary. However, in the presence of bending-extensional coupling, inplane behavior 
must also be considered, since it will occur even if the load is strictly transverse. According 
to [5], four different constraint cases involving various combinations of in-plane and out-of- 
plane displacement constraints and force resultant constraints can imply a simple support 
boundary condition; the condition used for the linear test problems solved herein is as 
follows: 

w = M n = N n - N nt = 0 

where n and t are used to designate directions normal and tangent to a boundary, respec- 
tively. For a formulation which excludes transverse shear deformation, the only condition 
to be imposed is the w = 0 requirement, since the constraints on moment and force resul- 
tants are considered natural boundary conditions and are therefore not enforceable in a 
Rayleigh- Ritz formulation. For a formulation which includes transverse shear deformation, 
the rotational quantity <j> n is also constrained for a simply supported edge; in the notation 
of the present work, for an edge where X\ = constant, (j> 2 = 0, and conversely for an edge 
where x 2 — constant, (f>\ — 0. 
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. Linear Polynomial Functions 



46 


Note: The terms “even” and “odd” refer to symmetry states as defined in Table 2.1. 
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The following two layup configurations can be combined with the two support/loading 
conditions to obtain four test problems: 

1. Orthotropic [0 4 /904] s , layer thickness = .005 in., size 20 in. by 9 in. 

2. Anisotropic with bending-extensional coupling: 

[±50/35]5, layer thickness = .00520087 in., size 19.25 in. by 6.5 in. 

For all four problems, comparisons were made with an in-house finite element code 
which included the effects of transverse shear deformation and employed mixed formulation 
finite elements [22]. For the orthotropic panels, a 12 by 8 grid of MD9-4 elements was used; 
these elements employ biquadratic shape functions in the approximation of the generalized 
displacements, and bilinear shape functions in the approximation of the stress resultants. 
For the anisotropic panels, a 12 by 6 grid of MD16-9 elements was used; these elements 
employ bicubic shape functions in the approximation of the generalized displacements, 
and biquadratic shape functions in the approximation of the stress resultants. The finite 
element model used is shown in Figure 4.2. 

To verify the correctness of the Rayleigh-Ritz solution, total strain energy and max- 
imum displacement were compared with the mixed formulation finite element solution, 
which is taken as the standard for comparison. Also, detailed examinations of the nonzero 
displacements were made by plotting each displacement at several points along the panel’s 
centerlines. 

4.2.1 Transversely Loaded Orthotropic Panel 

The transversely loaded orthotropic panel was used because it permits several simplifi- 
cations which greatly reduce problem size and complexity. First, approximation functions 
can be chosen which satisfy the boundary conditions; therefore, no penalty function is 
needed to enforce them. Second, solution of such a problem involves the use of a set of 
even symmetry state functions only; no odd symmetry state functions are needed, since the 
load and boundary conditions axe symmetric. Also, since there is no bending-extensional 
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coupling, no functions are needed to represent u i and u 2 . Taking these simplifications into 
account, the functions chosen to represent each displacement are: 


ui : 0 

u 2 : 0 



(TYl'KX\ ' 
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/"nnx2 \ 

cos 

( a , 

1 sin 

\ b ) 

• i 
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(Tl'KX'i \ 

sin | 

{ a J 

1 cos 1 
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• i 
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sin | 

V a ) 

1 sin ( 

< b ) 


(4.1) 


m, n = 1, 3, 5... 
m, n = 1, 3, 5... 
m, n = 1, 3, 5... 

Energy and center transverse displacement results are shown in Table 4.2. The column 
reporting number of terms refers to the number of terms representing each of the three 
displacements. Therefore, the total size of the linear system is three times the number of 
terms shown in Table 4.2. 


Table 4.2. Linear Results for [(L/QO^s, Transversely Loaded Panel 


Number 
Of Terms 

Energy 

w c , in. 

Value 

% Error 

Value 

% Error 

1 



4.799xl0~ 2 

2.97 

5 



4.651 xlO" 2 

0.22 

9 

0.1771 

0 

4.662xl0 -2 

0.03 

16 

0.1771 

0 

4.660xl0~ 2 

0.02 

25 

0.1771 

0 

4.661xl0“ 2 

0 

FEM 

0.1771 

— 

4.661xl0 -2 

— 


From the displacements shown in Table 4.2, the solution is essentially converged when 
five series terms are used. Examination of the plots shown in Figure 4.3 proves this 
assessment to be correct, since four of the five solutions are indistinguishable. Strain 
energy, however, reflects a slightly slower convergence, as well as a slight oscillation about 
the finite element solution as more terms are used. This behavior has proved typical of 
modified Rayleigh-Ritz technique. 


48 













4.2.2 Transversely Loaded Anisotropic Panel With 
Bending-Extensional Coupling 

The transversely loaded anisotropic panel was used because, unlike the orthotropic 
panel, it admits no simplifications, except that for a small enough load, it is a linear prob- 
lem. The presence of bending-extensional coupling means that approximation functions 
must be included for both «i and u 2 . Also, for the types of functions used, boundary 
conditions must be enforced using the penalty function method. Furthermore, both even 
and odd symmetry state functions must be included, since coupling between the two states 
is present. 

The full set of functions shown in Table 4.1 is used to approximate the response 
in this problem. Each displacement is approximated by a set of functions composed of 
the same number of even and odd symmetry state terms. For example, if 102 terms are 
used to approximate Uj, then the 51 even state functions consist of one linear term, 25 
sin( » 2 i£) co s( 1 ^) terms, and 25 cos(- 2J ^)sin(- 2 jp) terms, while the 51 odd state functions 
are of similar composition. Note that, according to Table 4.1. a, the functions representing 
w are an exception. In that case, since there is no even state linear term, only a total 
of 101 functions would be used. (If the even state linear term were included, then the 
cos(0) cos(0) term would have to be omitted, since it is also linear. The presence of both 
terms would introduce redundancy into the system, and the stiffness matrix would become 
singular and therefore unfavorable.) Also, because both inplane degrees of freedom are 
completely free, the odd state linear term representing u 2 is omitted. 

Energy and center transverse displacement results for the transversely loaded aniso- 
tropic panel are shown in Table 4.3. The column reporting number of terms refers to 
the number of terms representing each of the five displacements, with w and u 2 being 
represented by one less function, as discussed above. Therefore, the total size of the linear 
system is five times the number of terms shown in Table 4.3, minus two. 
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Table 4.3. Linear Results for [±50/35]5, Transversely Loaded Panel 


Number 
Of Terms 

Energy 

w c , in. 

Value 

% Error 

Value 

% Error 

18 

2.072xl0 -2 

0.99 

7.325xl0~ 3 

1.16 

38 

2.083xl0" 2 

0.46 

7.407xl0 -3 

0.05 

66 

2.087xl0" 2 

0.28 

7.403xl0 -3 


102 

2.089xl0 -2 

0.19 

7.405xl0 -3 


FEM 

2.093xl0~ 2 

— 

7.411xl0 -3 

E I 


From Table 4.3, the energy and transverse displacement are converged to three digits 
when 66 terms are used. Examination of the plots shown in Figure 4.4 confirms that using 
66 functions per displacement does indeed provide an acceptable solution, in that the 66 
term solution is nearly indistinguishable from the 102 term solution and from the finite 
element solution for the three larger displacements. Also, notice that good agreement is also 
obtained using the 38 term representation. (The finite element solution displays erratic 
results for displacements u\ and because they are much smaller than the maximum 
displacement in the panel.) 

4.2.3 Axially Loaded Orthotropic Panel 

The axially loaded orthotropic panel is another problem which permits several sim- 
plifications. First, since there is no out-of-plane behavior when the load is purely axial 
(and buckling has not occurred) no functions are needed to approximate w. Also, since the 
loading is symmetric, no odd state functions are needed. Therefore, the function set to be 
used consists of all the even state functions shown in Table 4.1, except that no functions 
are included to represent w or either of the two rotations <f>\ and <j> 2 . 

Boundary conditions for this problem are slightly more complicated than those for 
the transversely loaded panels. Because this panel, like all the rest, will exhibit coupling 
between axial and lateral behavior (i.e., Poisson effect), uniform displacement will not 
occur along the loaded edges. To prevent bowing of the loaded edges, the first derivative 
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of ui along those edges must be constrained to zero; therefore, the additional penalty term, 
consistent with the conventions of Section 3.2, is 


8h 


-I 


9K_ 

dx 2 


( 0 ,* 2 ) 


dJZ 1 

dx 2 


1 ( 0 , 12 ) 


X l2 dz 2 8X tl 


(4.2) 


Table 4.4. Linear Results for [CL/QO^s, Axially Loaded Panel 


Number 
Of Terms 

Energy 

8, in. 

u 2e ,in. 

Value 

% Error 

Value 

% Error 

Value 

% Error 

9 


0.02 

6.170xl0~ 4 

0.02 

1.433xl0 -5 

13.1 

19 


0 

6.170xl0~ 4 

0.01 

1.624xl0~ 5 

3.9 

33 

0.2777 

0 

6.171xl0" 4 

0 

1.501xl0 -5 

4.0 

51 

0.2777 

0 

6.171X10" 4 

0 

1.582xl0 -5 

1.2 

FEM 

0.2777 

— 

6.171 xlO" 4 

— 

1.564xl0 -5 

— 


Values of energy, loaded edge axial displacement 8 , and lateral displacement U 2 r along 
the simply supported edge are compiled in Table 4.4. The column reporting number of 
terms refers to the number of terms representing each of the two displacements; therefore, 
the total size of the linear system is twice the number of terms shown in Table 4.4. 

From Table 4.4, both the energy and loaded edge displacement are converged to three 
digits when only 9 terms are used. However, U 2 e never really converges, even when the 
full 51 terms are used. This behavior is due to the clamped boundary conditions along 
the loaded edge. Since ui is free on all four edges of the panel, a simple linear distribution 
is obtained, as shown in Figure 4. 5. a. However, since u 2 is constrained along the loaded 
edge, the Poisson effect which would normally cause a lateral expansion along that edge 
is inhibited. Therefore, since u 2 is not constrained along the unloaded edges of the panel, 
a look at the short centerline of the panel reveals a slightly more complicated nonlinear 
distribution for u 2 , as shown in Figure 4.5. b. 
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4.2.4 Axially Loaded Anisotropic Panel With 
Pending- Extensional Coupling 

The function set used for this problem is identical to that used for the transversely 
loaded orthotropic panel, except that the odd state linear term is included in the repre- 
sentation of U 2 ■ The “straight edge boundary condition” introduced for the axially loaded 
orthotropic panel is also used here, and prevents not only bowing due to the Poisson effect, 
but also edge shearing due to shear-extensional coupling, as shown in Figure 2.2. 

Results for the axially loaded anisotropic panel problem are shown in Table 4.5. Since 
the full set of functions shown in Table 4.1 must be used to approximate the response, the 
number of terms reported in Table 4.5 must be multiplied by five, and then one subtracted, 
to obtain the size of the full linear system. 


Table 4.5. Linear Results for [±50/35]5, Axially Loaded Panel 


Number 
Of Terms 

Energy 

8 , in. 

Value 

% Error 

Value 

% Error 

18 

wmgm 

5.41 

1.684xl0“ 3 

5.34 

38 

II 

3.01 

1.727xl0" 3 

2.92 

66 

WSM 

2.04 

1.745xl0 -3 

1.91 

102 

HI 

1.51 

1.754xl0~ 3 

1.41 

FEM 

■H 

— 

1.779xl0 -3 

— 


From Table 4.5, it is evident that this problem is the most demanding of the four, in 
terms of the number of functions needed to obtain an accurate solution. Three-digit accu- 
racy is never obtained in any of the solutions shown. Examination of Figure 4.6 indicates 
that a 102 term approximation provides a good estimate of the finite element solution, 
while a 66 term approximation gives results very nearly as good. Note that in-plane be- 
havior is predicted well using only a few terms, while the out-of-plane quantities are more 
difficult to model. (The finite element solution displays erratic results for displacements 
<f>i and <f > 2 because they are much smaller than the maximum displacement in the panel.) 
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4.3 Effect of the Penalty Parameter on Accuracy of Linear Results 

For the results just presented, no mention was made of the penalty parameter; how- 
ever, choice of this parameter can have a profound effect on the accuracy of the solutions. 
Too large a penalty parameter can result in an ill-conditioned equation system, while a 
choice at the other extreme can mean that boundary conditions are not adequately satis- 
fied. In this section, results of parametric studies of the penalty parameter axe presented 
to give the reader some appreciation of the effect of this parameter. In this chapter, an 
overall penalty parameter of 10 8 was used, unless otherwise noted. 

First, examine Table 4.6. This table presents results from solution of the axially loaded 
anisotropic panel of Section 4.2.4 for different values of the penalty parameter r. Several 
different parameters are presented; the following list defines each one: 
r = Penalty parameter 
U = Strain energy 

\R\ — Magnitude of the residual vector 12, calculated by substituting the Ritz 
coefficients back into the linear equilibrium equation KijX.j — Pi. 

8 = End shortening at the center of the panel’s loaded edges 
8 BC = Number of digits to which the loaded edges are straight 
w c = Center transverse displacement 

w BC = Order of magnitude satisfaction of the w = 0 constraints along the 
panel’s four edges 

U 2 BC = Order of magnitude satisfaction of the U 2 = 0 constraints along the 
panel’s four edges 

It is evident from these results that all three of the response variables, 27, 8, and 
w c , are essentially insensitive to variation of the penalty parameter when the penalty 
parameter is at least 10 8 , and the residual is at an acceptable level regardless of the value 
of r. However, boundary condition satisfaction is much more sensitive to variations in the 
penalty parameter. 

It seems that the easiest boundary condition to satisfy is the w = 0 condition; this 
condition is at an acceptable level for any value of r greater than 10 7 . The most difficult 
condition to satisfy appears to be the U 2 = 0 condition along the loaded edge; this condition 
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Table 4.6. Penalty Results - Overall Penalty Parameter* 


r 

u 

1*1 

8xl0~ 2 

8 BC 

C4 

! 

O 

r-H 

X 

0 

3 

w BC 

u 2 BC 

(u!? ax ~ 10~ 4 ) 

10 5 

.6011 

10" 15 

1 

1 

-.1131 

10~ 6 

10- 4 

10 6 

.5810 

10" 15 


1 


10- 7 

10~ 5 

10 7 

.5707 

10-14 

.1758 

2 


10~ 8 

10" 5 

10 8 

.5670 

10- 13 

.1745 

3 

-.1053 

10~ 9 

10" 6 

10 9 

.5662 

10-12 

.1742 

4 

-.1052 

10~ 10 

10- 7 

10 10 

.5661 

10- 11 

.1742 

5 

-.1052 

10' 11 

10~ 8 

10 u 

.5661 

IQ- 10 

.1742 

6 

-.1052 

10~ 12 

10- 9 


* See text for definitions of variables presented in this table 


does not become sufficiently small until r reaches 10 11 . (Sufficiently small is defined as 
six orders of magnitude smaller than the largest displacement in the panel, which is 8, at 
10 _3 .) Furthermore, it is also difficult to obtain a constant 8 condition along the loaded 
edge; this is not accomplished to six digits until r reaches 10 11 . 

More insight into these difficulties are gained by examining the contour plots shown 
in Figure 4.7. Look first at the contour of ui; from this plot, it is evident that the effect 
of an unsatisfied straight edge condition is that the edge is slightly bowed. Note that the 
skewing caused by the presence of shear-extension coupling is successfully prevented along 
the loaded edge; however, one cannot blame this coupling for the difficulty in satisfying 
the straight edge condition, since similar problems were encountered in working with the 
[±30/90] 5 panel, which has no extension-shear coupling. 

Rather, as is evident from the u 2 contour plot, the cause of the difficulty is the 
requirement that u 2 be zero along the loaded edge. Notice that steep u 2 gradients are 
present in the corners of the panel; since u\ and u 2 are very strongly coupled via the 
well-known Poisson effect, there is a strong tendency for the loaded edge to bow inward, 
and for the corners to also move inward. Therefore, much higher penalty parameters are 
required in order to satisfy both the constant 8 and u 2 = 0 conditions. 
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In order to put these difficulties in perspective, it is important to note that, even 
though the boundary conditions are not satisfied to as many digits as might be desirable, 
the desired character of the response is obtained. That is, the gradients in 112 at the 
corners are in fact present, as is evident from the contour plots, and skewing caused by 
shear-extensional coupling is totally prevented along the loaded edges. 

In addition, experience with nonlinear analysis of these problems indicates that it 
is unwise to impose overly stringent precision requirements on boundary conditions. For 
example, orthotropic problems were solved in which an overall penalty parameter of 10 9 was 
used, with a penalty parameter of 10 11 being used to enforce the straight edge boundary 
condition. When this problem was run on a VAX 11/785, which carries sixteen significant 
digits in all calculations (in double precision), acceptable results were obtained. However, 
the same problem with the same penalty parameter was then solved on a Cray-2, which 
carries only fourteen significant figures (in single precision); in that case, convergence was 
not obtained even on the first load step, due to the ill-conditioning cause by the large 
penalty parameters. Conversely, when the penalty parameter for the entire panel was 
decreased to 10®, converged solutions were easily obtained on both machines. Furthermore, 
note that since only half as many functions were required to solve this orthotropic problem 
than are needed for an anisotropic problem, ill-conditioning due to the size of the penalty 
parameter can only get worse for anisotropic problems. Therefore, even though Table 4.6 
indicates that acceptable linear solutions can be obtained using high penalty parameters, 
care must be exercised when choosing penalty parameters for nonlinear problems. 

4.4 Linear Results Using Reduced Basis Techniques 

From the results of the previous sections, it is evident that acceptable solutions can 
be obtained using subsets of the functions shown in Table 4.1. However, as indicated by 
results for the axially loaded anisotropic panel, the system of linear equations required for 
an accurate solution can become quite large; in that case, 509 equations were required. The 
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reduced basis technique described in Section 3.4 is one way to reduce the size of the equation 
system. Recall that this technique involves approximating the solution using a set of basis 
vectors which were essentially sensitivity derivatives of the response with respect to the 
anisotropic material coefficients, evaluated at the point where the anisotropic coefficients 
are zero. 

In this section, solutions obtained using this reduced basis technique will be examined 
for the axially loaded anisotropic panels of Section 4.2.4. In particular, two details of 
these results are of interest. First, it is of course important to determine how well the 
chosen basis vectors approximate the full system. This determination is easily made by 
examining local and global quantities for the panel being analyzed. Secondly, the basis 
vectors themselves should be examined in order to determine whether they are linearly 
independent. 

In order to measure the linear independence of a given set of vectors, the Grammian 
may be used. As explained in [23], the Grammian is the determinant of the Gram matrix, 
which is formed as follows: 

M cd = T cI T Id (4.3) 

That is, elements of the Gram matrix are formed by taking appropriate dot products 
involving the basis vectors. If the basis vectors are normalized beforehand so that the 
length of each vector is one, then the Grammian will range between zero and one. If the 
Grammian vanishes, then the basis vectors can be shown to be linearly dependent, while 
a Grammian of one indicates complete linear independence. Note that, as will be shown 
herein, the Grammian can be used to determine the required number of basis vectors, thus 
eliminating the need for detailed studies of response quantities. 
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4.4.1 Accuracy in Approximating the Full System 


Shown in Table 4.7 are several key quantities calculated using reduced systems ranging 
in size from one to fifteen basis vectors. These results axe for the axially loaded anisotropic 
plate, whose full model involves 329 equations. From this table, both the strain energy and 
the end shortening of the panel converge very rapidly; only four basis vectors axe required 
for accurate representation of both quantities. However, in order to obtain a negligible 
residual vector, at least eight basis vectors are required. The effect of this last result is 
seen when smaller displacements, such as the maximum transverse displacement w c , are 
calculated; at least seven basis vectors axe required before w c is converged to six digits. 
Therefore, accurate solution of this problem requires the use of at least eight basis vectors. 


Table 4.7. Variation of Results with Number of Basis Vectors Used 


X 

6 

U 

w c 

i*i 

1 

.161227xl0“ 2 

0.523919 

-9.29727xl0 -4 

10° 

2 

.174011 xl0~ 2 

0.565462 

-1.00345xl0 -3 

1 

O 

H 

3 

. 1 74465 x 10 ~ 2 

0.566938 

-1.04855xl0~ 3 

10~ 2 

4 

.174469 xl0~ 2 

0.566950 

-1.04893xl0 -3 

10~ 2 

5 

(converged) 

(converged) 

-1.05267xl0 -3 

10~ 3 

6 

. 

. 

-1.05276xl0 -3 

1 

o 

i-H 

7 

- 

- 

-1.05268xl0 -3 

I- 4 
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1 

0» 
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(converged) 

10~ 6 
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. 

10- 7 
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1 
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12 
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4.4.2 Grammian Results 

Table 4.8 lists values of the Grammian for different numbers of basis vectors. Again, 
these results are for the anisotropic axially loaded panel, and 329 functions were used in 
the full model. These results are consistent with those of the previous subsection since 
calculation of more than eight basis vectors causes the Grammian to essentially vanish. 
Furthermore, the results of Table 4.8 indicate that, for this class of problems, the Grammian 
provides a good prediction of the number of basis vectors required, and therefore alleviates 
the need to examine detailed results in order to select an appropriate reduced system size. 














Other numerical experiments have shown that, as one might expect, when fewer ap- 
proximation functions are used and the number of unknowns is therefore smaller, fewer 
basis vectors are needed to approximate the reduced number of unknowns. This is because 
the space defined by the smaller set of unknowns is of lower dimension, and so fewer basis 
vectors are needed to span it. Shown in Table 4.8.b are values of the Grammian calculated 
for the 89 function model; in that case, only seven basis vectors are needed. 

4.5 Nonlinear Results 

The problems for which nonlinear results are presented involve axially loaded or- 
thotropic and anisotropic panels similar to those presented in the previous sections, as 
well as a 6.5 inch wide, axially loaded [±30/90 ] 5 panel similar to that described in Chap- 
ter 5. Note that the dimensions of the orthotropic panel are different from those of the 
anisotropic panels; this was done in order to obtain a control specimen which, like the 
two anisotropic panels, buckled into a w-symmetric configuration with respect to the short 
centerline of the panel. 

The two nonorthotropic panels are identical to the corresponding panels discussed in 
Chapter 5 in both lamination and thickness; however, their widths and lengths have been 
modified to compensate for the slight difference in the way boundary conditions must be 
modeled in the Ritz and finite element analyses. In the finite element models of Chapter 5, 
an extra ring of elements was added to simulate the finite-width supports used to impose 
boundary conditions in the experiments. However, in the Ritz analysis, the supports are 
modeled by simply shortening both the panels’ dimensions by amounts corresponding to 
the extra ring of elements added in the finite element model. All five degrees of freedom are 
constrained along the loaded edges, while only w is constrained along the other two edges. 
(In the linear problems presented earlier, both w and <f>\ were constrained along the long 
panel edges.) Note that there is no reason that the exterior penalty function approach 
presented earlier could not be used to enforce the more complex boundary conditions; 
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however, it was decided that for purposes of this work, the additional complexity was 
unnecessary. As a result, quantitative comparison is not made between the experimental 
results and the other results presented in this section; however, qualitative results are 

I 

discussed. 

Also, in Chapter 5, load was introduced at one end of the panel while the other end 
remained motionless; for the problems of the present section, load is introduced at both 
ends of the panel. The two loadings are functionally identical; the latter was used in 
order to obtain a symmetric loading condition, which is required by the Rayleigh-Ritz 
formulation. Schematics of the Rayleigh-Ritz analysis models are presented in Figure 4.1. 

Three other sets of results are presented for comparison with the Ritz results; these are: 
STAGS results, results from the mixed formulation finite element model used previously 
for comparison with the linear results, and the experimental results presented in Chapter 
5. The STAGS results presented in the following sections were obtained using models 
which correspond in geometry, boundary conditions, and lamination with the Ritz models, 
and so will be slightly different from the results presented in Chapter 5. Results obtained 
from the mixed formulation finite element analysis also correspond to the Ritz models 
in geometry, lamination and boundary conditions. (For all three analyses, the mixed 
formulation analyses employed the MD9-4 finite element described in Section 4.2.) Also, 
recall that both the Rayleigh-Ritz results and the mixed formulation results include the 
effects of transverse shear deformation, while the STAGS results do not. 

The solution strategy chosen for the Rayleigh-Ritz analysis of the nonlinear post- 
buckling problem is the well-known Newton-Raphson technique. This technique, although 
computationally expensive compared with most transformation-type methods, is simple 
to implement and gives accurate results. The particular solution procedure used herein 
consists of the following steps: 

1. Obtain an approximation for the buckling load of the panel being analyzed. This was 
accomplished using the finite element code STAGSC-1, since no eigenvalue capability 
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was included in the Rayleigh-Ritz code. (See Chapter 5 for further discussion of 

STAGSC-1.) 

2. Solve the linear problem at an initial load of 70% to 90% of the buckling load; this 
solution is an initial estimate for the first nonlinear load step. 

3. Form the nonlinear Newton- Raphson equations as described in Chapter 3; the solution 
from the previous iteration (or the linear problem if this is the first load step) is used 
in the formation of these equations. 

4. Solve this linear set of equations for a correction AX; to the previous solution X;_!. 

5. Form the new solution vector X; = X;_i + AX;. 

6. Calculate an error measure 

1 [AX;| 
e n X; 

where n is the total number of degrees of freedom. When this error is less than 10 ~ 6 , 
the solution is considered converged. If the solution is not converged, return to step 
3. 

7. When a converged solution has been obtained at the current load step, the load may 
be incremented. The size of this incrementation depends on the proximity of the 
current load to the buckling load. If the load is within five to ten percent of P cr , then 
incrementation proceeds very slowly. At other points in the solution procedure, the 
load may be incremented at larger intervals. 

Throughout this procedure, an imperfection must be imposed in order to trigger buckling; 
without it, no branching to the stable secondary path will occur. The imperfection is 
achieved by applying a sinusoidal transverse load with the same number of half-sine waves 
in both directions as the first buckling mode shape. (It happens that in all the cases 
presented, the desired response corresponds to the first buckling mode.) In the case of the 
STAGS analyses, imperfections of exactly the same shape as the buckling mode shape were 
used, just as in Chapter 5, and in the case of the mixed formulation results, no imperfection 
was needed. Therefore, considerable differences between the three sets of results will be 
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seen in the immediate neighborhood of buckling; however, imperfection related differences 
will disappear as the load or displacement is increased beyond buckling. For a more general 
discussion of the use of imperfections to trigger buckling, as well as the entire nonlinear 
solution procedure, refer to Chapter 5. 

For the Rayleigh- Ritz analyses, the same sets of functions are used to solve the post- 
buckling problems as were used for the linear problems, with only one exception: all five 
degrees of freedom (u\, 112 , <f> 1 , <f> 2 , and w) must be represented for all problems, regard- 
less of the laminate, since both inplane and out-of-plane deformations occur. Results are 
presented for two problem sizes: 38 and 66 terms per displacement degree of freedom, 
resulting in full systems of 189 and 329 equations, respectively (since w is represented 
by one less function than are the other four displacements). Larger problems would have 
been desirable, but problem size was limited by the computer hardware. In the case of the 
orthotropic problem, half as many functions are used, since only the even function set is 
required, resulting in full systems of 94 and 164 equations. 

In the case of the STAGS analyses, the same 411-type elements that were used for 
the analyses of Chapter 5 were used here. In all three problems, a uniform grid was 
used, with the loading and boundary conditions applied as described in Chapter 5. For 
the orthotropic problem, a discretization resulting in 10 by 12 elements per half-wave of 
the buckling mode shape was used, for a total of 360 elements, while for the orthotropic 
problems, the discretization gave 10 by 10 elements per half-wave, resulting in a total of 
300 elements. Studies were made to assure that these levels of discretization provided 
converged solutions. 

Data presented in the plots of the following sections has all been normalized with 
respect to the panel buckling load. In the case of the STAGS, Rayleigh-Ritz, and mixed 
formulation finite element results, this means normalization with respect to the buckling 
loads obtained via a STAGS bifurcation buckling analysis. The experimental results shown 
were normalized with respect to the buckling loads presented in Table 5.6. 
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It is important to note that all the Ritz solutions presented in the following sections 
were solved on a Cray-2 supercomputer. Because Cray-2 single precision is slightly less 
accurate than than the double precision available on a VAX 11/785, careful choice of an 
appropriate penalty parameter becomes more important, as was noted earlier. In all three 
problems, the penalty parameter was set to 10 8 . 

4.5.1 Orthotropic Panel 

The orthotropic panel analysis was taken to a load equal to about three times the buck- 
ling load, as determined by bifurcation buckling analyses done with STAGSC-1. Shown 
in Figure 4.8 is a plot of load versus end- shortening for this panel which includes four of 
the five sets of results described in the previous section, since no experimental results were 
available for the orthotropic panel. (For reference, the buckling load for this panel was 
1452 lb.) 

Aside from the expected differences between the results in the immediate neighborhood 
of buckling, the primary differences between the three sets of results are seen late in 
postbuckling. The STAGS results are slightly stiffer than the other three sets of results, 
which agree well during the early portion of postbuckling. Since the 33 term Ritz solution 
agrees well with results from the mixed formulation finite element analysis, that solution is 
considered converged. The 19 term solution, while not in total agreement with the mixed 
formulation finite element solution, still provides an accurate solution early in postbuckling. 

Since the only real difference between the STAGS and Ritz analysis is the presence of 
transverse shear flexibility in the Ritz analysis, this difference in stiffness late in postbuck- 
ling could be attributed to transverse shear deformation. Since agreement is good between 
the mixed formulation results and the Ritz analysis, the claim that transverse shear defor- 
mation is the cause of the disparities in stiffness appears to be well supported. However, 
note that the degradation in stiffness, as defined in Chapter 5, is essentially identical for 
the two solutions, due to slight differences in initial postbuckling stiffness. 
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The agreement between the 33 term Ritz solution and the mixed formulation finite 
element solution is also seen on a pointwise basis, as shown in Figure 4.9, a plot of center 
transverse displacement. However, the 19 term Ritz solution is evidently not a good 
estimate of the behavior, and similarly for the STAGS solution. The extreme disagreement 
between STAGS and the other results in the immediate neighborhood of buckling is due 
to the larger imperfection imposed in the STAGS analysis, as noted earlier. 

Finally, examining the two contour plots shown in Figure 4.10 (one from the Ritz 
analysis and one from the STAGS analysis), it is also evident that overall qualitative 
agreement between the two analyses is good for this problem. 

4.5.2 [±30/90]* Panel 

The analysis of the [±30/90] 5 panel was also taken to a load equal to about three 
times the buckling load. Shown in Figure 4.11 is a plot of load versus end-shortening for 
this panel which includes all five sets of results described earlier. From this plot, it is 
evident that for the Ritz analysis, 38 functions per displacement degree of freedom is an 
inadequate representation of the problem; however, agreement is considerably better for 
the 66 function case. Even so, this higher level of approximation is still not sufficient, 
since agreement with the mixed formulation solution is not attained. (The fact that the 
Ritz solution agrees well with the experiment is not taken as proof of convergence, since 
boundary conditions were not modeled precisely.) 

The plot of center transverse displacement, shown in Figure 4.12, displays a very dif- 
ferent character from the end shortening results. Apparently, pointwise agreement between 
the different solutions is not attained for this panel. The analysis which comes closest to 
the experiment is the STAGS analysis; then, being next in order of agreement, both the 
mixed formulation results and the Ritz analysis are less stiff than both the experiment and 
the STAGS results. One might suspect that the large imperfection used in the Ritz anal- 
ysis is the cause of the disparity for that set of results; however, runs were made to ensure 
that the postbuckling solution did not vary significantly when the imperfection was made 
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an order of magnitude larger, thus eliminating imperfection from consideration. Rather, 
the most likely cause of the disparity between the mixed formulation results and the Ritz 
results is simply a lack of convergence on the part of the Ritz analysis. Further discussion 
of the differences between the experimental results and the STAGS analysis may be found 
in Chapter 5. 

Finally, examining the two contour plots shown in Figure 4.13 (one from the Ritz anal- 
ysis and one from the STAGS analysis), qualitative agreement between the two analyses is 
attained for this problem in that both analyses display a three half-sine-wave postbuckled 
shape. However, the Ritz analysis, presumably due to its lack of convergence, displays 
some distortion near the center of each half wave. 

4.5.3 [±50/351. Panel 

The analysis of the [±50/35] 5 panel was taken to a load equal to about two and one- 
half times the buckling load. Shown in Figure 4.14 is a plot of load versus end- shortening 
for this panel which again includes all five sets of results described earlier. From this plot, 
it is again evident that for the Ritz analysis, 38 functions per displacement degree of free- 
dom is an inadequate representation of the problem, while the 66 function approximation 
exhibits better agreement with the other solutions, but still is not adequate. Excluding 
the Ritz solutions and examining the three remaining solutions late in postbuckling, the 
stiffness degradation seen in the experiment is not predicted by either the STAGS analysis 
or the mixed formulation results. Therefore, one cannot conclude that transverse shear 
deformation is responsible for postbuckling stiffness degradation in this problem. 

The plot of center transverse displacement, shown in Figure 4.15, displays a character 
much different from that seen for the [±30/90]s panel, in that better agreement is obtained 
for all results except the 38 term Ritz solution. Obviously, the 38 term approximation is 
totally inadequate. The disparity between the 66 term Ritz solution and the other curves 
again indicates the lack of convergence in the Ritz solution. 
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Finally, examining contour plots for this panel, shown in Figure 4.16, qualitative 
agreement between the Ritz and STAGS analyses is attained for this problem, since both 
solutions exhibit a three half-sine-wave postbuckled shape, as well as the skewing of the 
entire pattern which is characteristic of panels with shear-extensional coupling. As with 
the [±30/90 ] 5 panel, the Ritz solution displays some distortion near the centers of the 
three half-waves. This is particularly true for the center half-wave, where the oblong shape 
predicted by the STAGS analysis is rotated ninety degrees in the Ritz analysis. 
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Chapter 5 

Experimental and Numerical Studies Of 
Buckling and Postbuckling Behavior of Panels 

5.1 Characteristics of Test Specimens 

Two of the primary aims of the present study axe: 

• To investigate the buckling and postbuckling responses of rectangular laminated com- 
posite panels possessing anisotropy and bending-extensional coupling 

• To use readily available finite element modeling techniques for the analysis of such 
panels 

In order to collect experimental data which quantifies the response of realistic lam- 
inates possessing both anisotropy and bending-extensional coupling, six graphite-epoxy 
panels were constructed; three of the panels were made using a [±50/35]s layup config- 
uration, and the other three were constructed using a [±30/90]5 layup. The [±50/35 ] 5 
laminate has been used in the design of a prototype aircraft known as the HiMAT (Highly 
Maneuverable Aircraft Technology) vehicle, while the [±30/90]s laminate was chosen as 
a control specimen, since it has no membrane anisotropy. Note that both of these layups 
are unsymmetric and have a total of fifteen layers. The two sets of panels exhibit all the 
various couplings described in Chapter 2 to differing degrees. 

The first layup, [±50/35]s, possesses all the different couplings; its constitutive matrix, 
which corresponds exactly to the matrix C defined in Chapter 2, is shown in Table 5. 2. a. 
The second layup, [±30/90]5, exhibits no extension-shear or transverse shear coupling, 
a lower degree of bending-twisting coupling, and a higher degree of bending-extensional 
coupling. The constitutive matrix for this laminate is shown in Table 5.2.b. All the panels 
were made from Hercules AS4-3502 graphite-epoxy unidirectional preimpregnated tape. 
Properties of this material system are listed in Table 5.1. (Note that the matrices shown 
in Table 5.2 were calculated using a panel thickness of 0.08 inches. The actual panels 
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Table 5.1. Material Properties for Hercules AS4-3502 Graphite/Epoxy Tape 


E c n 

18.5xl0 6 

77»C 

1.64xl0 6 

G\2 

0.832xl0 6 

Vl2 

0.35 


had thicknesses slightly different from this nominal value, and so the constitutive matrices 
used in the analyses were slightly different from those shown in the table; the purpose of 
presenting the nominal matrices is to show the approximate magnitude of the coupling 
involved.) 

One other parameter, the aspect ratio, was also varied. Three panels in each of the 
two layup sets were cut to widths of 4, 5.5, and 7 inches, with all lengths being a constant 
20 inches. Average thicknesses were calculated for each panel based on eight measurements 
taken over each panel’s surface. These average thicknesses have been used in all analyses 
presented herein. Widths and thicknesses of each panel are summarized in Table 5.3. 

The six panels were all tested to failure, and each test was modeled using the finite 
element code STAGSC-1 (hereafter referred to as STAGS and described fully in reference 
[24]). Among other features, STAGS possesses a well-tested set of nonlinear solution 
algorithms, as well as a linear bifurcation buckling analysis capability, both of which were 
useful in the present study. The plate bending elements included in STAGS do not permit 
transverse shear deformation. Correlation between experimental and finite element results 
is examined in detail. 

5.2 Description of the Experiment 

5.2.1 Test Setup. Measurement Techniques, and Procedure 

A schematic and photograph of the overall test setup is shown in Figure 5.1; this 
setup is essentially the same for all the tests conducted in the present study. Tests were 
conducted using a 120 kip hydraulic test machine. The loading was axial compression in 
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Table 5. 2. a. Constitutive Matrix for [±50/35]s Panels, Thickness = 0.08" 


’ .480977 x 10* .322167 x 10* 

Cjj i, j 
.132166 x 10* 

= 1, 2, ... 6 
-.593526 x 10 3 

.430193 x 10 2 

-.219893 x 10 3 - 

.483230 x 10* 

.684693 x 10* 

.430193 x 10 2 

.507488 x 10 3 

.268056 x 10 3 


.342118 x 10* 

-.219893 x 10 3 

.268056 x 10 3 

.430193 x 10 2 

Sym. 

r 


.226447 x 10 3 

■ A » 4 1 

.150944 x 10 3 
.225668 x 10 3 

.652591 x 10 2 
.357176 x 10 2 
.160296 x 10 3 . 


,536847 x 10 s .272511 x 10 4 

.272511 x 10 4 .537153 x 10* 


Table 5.2.b. Constitutive Matrix for [±30/90] 5 Panels, Thickness = 0.08" 


C ij i, j = 1, 2, ... 6 


r .613178 x 10* 


.191093 x 10* 
.613178 x 10* 


0 . 

0 . 

.211043 x 10 6 


.122973 x 10 4 
.371607 x 10 3 
.676815 x 10 3 
.285378 x 10 3 


.371607 x 10 3 
-.197295 x 10 4 
.247720 x 10 3 
.889553 x 10 2 
.290716 x 10 3 


*. 3 = 7 , 8 


.676815 x 10 3 ‘ 
.247720 x 10 3 
.371607 x 10 3 
.338407 x 10 1 
.123860 x 10 1 
.983070 x 10 2 . 


Sym. 


Co = 


' 0.537000 x 10* 0. 

0. 0.537000 x 10* 


the long panel direction, and consists of an applied displacement at the bottom edge of 
the panel. 

The boundary conditions were imposed by steel supports attached to the panels along 
all edges. (See Figure 5.2.) The panels’ short edges were clamped along their entire length; 
each clamping fixture covered an area along the panels’ short edge extending from that 
edge to a line 0.375 inches into the interior of the panel, for a total clamped area of w 
by .375 inches, where w is the panel width. From a modeling standpoint, these clamping 
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Table 5.3. Measurements of Tested Panels 




Average 

Layup 

Width, in.f 

Thickness, inxlO -3 


4 

76.788 

[±30/90] 5 

5.5 

77.488 


7 

78.763 


4 

77.313 

[±50/35] 5 

5.5 

77.150 


7 

78.013 


f All panels are 20 inches long. 


supports imply, in the ideal case, that all displacements and rotations are constrained to 
zero in the area covered by the supports. 

The long edges of the panel were constrained using steel knife-edge supports, as shown 
in Figure 5.2. These supports extended only .25 inches into the interior of the panel. They 
constrained the out-of-plane displacement ( w ) to zero along their line of contact .25 inches 
inside the panel edge; however, rotations, as well as inplane displacements ui and « 2 , were 
free, since these supports did not contact the clamping supports. 

Five kinds of measurements were made during the tests: load magnitude, out-of-plane 
displacement, displacement of the loaded edge (or head), strain magnitudes, and qualitative 
characterization of the overall buckle pattern. Load was measured by the test machine 
load cell, and was displayed throughout the test on a properly calibrated voltmeter. Strain 
measurements were made using foil- resistance strain gages placed at various locations on 
the panel. Out-of-plane displacement was measured using induction- type spring/plunger 
transducers known as direct current differential transducers, or DCDT’s. Displacement of 
the loading head was also measured using a DCDT, and the overall buckle pattern was 
displayed using the moire fringe technique. Strain gage and DCDT voltages were recorded 
on magnetic tape using an automated data acquisition system and were translated to strain 
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units and inches, respectively, by postprocessing software. Load data was also recorded in 
this way. 

The strain gage/DCDT layout was decided on the basis of a preliminary STAGSC-1 
linear bifurcation-buckling finite element analysis. This analysis provides a good estimate 
of buckling load and a reliable general indication of the buckle patterns to be expected, and 
has proven to be accurate enough for final buckling load calculations. The level of accuracy 
exemplified by these buckling load predictions indicates that the prebuckling behavior for 
the panels studied was not appreciably nonlinear. Bifurcation analyses have also been 
helpful in indicating closely spaced modes so that panels may be gaged in case either the 
first or the second mode is observed. (Closely spaced modes occurred in three of the six 
panels tested.) 

Once the buckle pattern is known from the bifurcation analysis, gages are placed to 
record the maximum expected strains and displacements in the panel (e.g., in the middle 
of a half-wave) and also at the center of the panel, as a point of reference. (On two of the 
panels, strain gages were also placed on nodal lines to investigate suspected failure mech- 
anisms.) The same strain gage pattern is used for both sides of the panel, so that there 
are pairs of gages mounted back-to-back over the surface of the plate. This arrangement is 
useful since bending behavior is easily observed by the divergence of back-to-back strains. 
Divergence occurs because when the panel bends, one of the gages undergoes compression 
while the other experiences tension. (See Figure 5.3.) In the case of a perfect panel with 
little or no bending-extensional coupling, the point of divergence is coincident with buck- 
ling, whereas a panel with a high degree of bending-extensional coupling or imperfection 
exhibits divergence at loads far below buckling. 

In order to observe the buckle pattern as the test progresses, the moire fringe technique 
is used. This technique consists of first painting one side of the panel white so that it is 
highly reflective. A piece of transparent plastic film covered with very thin vertical lines (50 
per inch) is placed in front of the panel’s painted side. Then, a high-intensity light source 
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is directed at the panel throughout the test. The variation in out-of-plane displacement, 
over the panel surface during buckling and postbuckling makes it possible for the observer 
to see a contour plot of the panel’s out-of- plane displacement during the entire test. This 
technique is extremely helpful to the engineer both during and after the test in determining 
the overall behavior of the plate without taking reams of data. For example, the effect of 
anisotropy (in particular, extension/shear coupling) is very obvious, since one can actually 
see the buckle pattern skewing during the test. 

To study various failure characteristics of the six panels, each panel was tested twice; 
the only difference between the two separate runs was the stopping point of the test. The 
stopping point of the first test is determined by sounds the panel makes as it is loaded. 
For the material tested, impending failure is easily detected as a series of cracking or 
popping sounds; the first audible crack was then used as the stopping point for the first 
run. Although admittedly an imprecise technique, this method generally allows one to 
examine the area of failure initiation long before the panel is totally destroyed. 

The point of crack initiation is located using a nondestructive test technique known 
as the C-scan. The C-scan is an ultrasonic technique, in which the panel is submerged 
in water, sound waves are passed through the panel, and their attenuation is recorded at 
several thousand points over the panel’s surface. Maximum attenuation is seen as a white 
area on the C-scan, while minimum attenuation is indicated by black areas. At points 
where the panel has become delaminated, sound does not penetrate the panel, and so the 
initial failure is identified as a white area on the C-scan, while black areas indicate that 
the material is still structurally sound. (See Figure 5.4 for a sample C-scan; the remaining 
C-scans are collected in Appendix C.) Unfortunately, the C-scans also show the strain 
gages and their associated wiring, since the gages had to be left intact for the second test; 
however, failure initiation is still clearly visible in most of the C-scans. After C-scanning, 
the second test run was conducted, in which the panel was allowed to fail totally, so that 
the final failure location, load, and other data could be recorded. 
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5.2.2 Observations Regarding Experimental Techniques 


During this study, important differences between the test methods for orthotropic and 
anisotropic panels were noted. In tests involving orthotropic panels, two strain symmetry 
conditions can be used to ensure that the loading head is level. The first symmetry 
condition is that back-to-back gages must register equal strains at some low load prior 
to buckling (about 700 lbs. for these panels). The assumption is that if the loading 
head is level in the plane perpendicular to the plate, then the load is purely axial and 
there is no transverse component and therefore no bending. Secondly, gages in symmetric 
locations with respect to the horizontal centerline should show equal strain; this indicates 
that the loading head is level in the plane of the plate. (Refer to Figure 5.5.) For panels 
with bending-extensional coupling, the first requirement is never satisfied except at zero 
load, since bending is experienced whenever an axial load is applied. Furthermore, the 
symmetry conditions for anisotropic panels are not the same as those for isotropic panels. 
For anisotropic panels, inversion symmetry rather than mirror symmetry is observed, as 
described in Chapter 2. 

Before these differences were realized, very poor prebuckling stiffness agreement was 
obtained (30 to 40% error); however, once different methods were used to ensure a level 
loading head, error was reduced to 2% or less. To level the loading head, a very simple 
procedure is used. The loading head and platen are brought together and an extremely 
high load (50,000 - 60,000 pounds) is applied to the loading head. Then, all the leveling 
screws are tightened. After performing these simple steps, the loading head and platen axe 
parallel. 
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5.3 Finite Element Model 

5.3.1 Discretisation and Boundary Conditions 

All the experiments were modeled using the STAGSC-1 finite element code. Sample 
runstreams for both linear bifurcation buckling and nonlinear analyses are shown in Ap- 
pendix D. Figure 5.6 is a plot of the undeformed model geometry for the 5.5 inch wide 
plates, with boundary conditions as noted. In addition to the edge boundary conditions 
shown, additional constraints were placed on interior nodes to facilitate modeling of the 
supports described in Section 5.2.1. These constraints, depicted in Figure 5.6, are applied 
using STAGS load definition cards. Also, although the experiments are conducted using 
an applied displacement at one edge of the plate, the problem is modeled using an applied 
load at one point on the loaded edge, in conjunction with multi-point constraints to ensure 
that the loaded edge remains straight. This technique is easier than formally applying 
a displacement, since it makes postprocessing for equilibrium forces unnecessary. All the 
analyses were done using the STAGS 411 plate element. 

As can be seen from Figure 5.6, the grid is fairly fine and is nonuniform but mirror- 
symmetric about the panel centerlines. The nonuniformity arises from the need to compare 
STAGS displacements with experiment data at certain points on the panels. Symmetric 
grids were used because deviation from a symmetric pattern adversely affects results unless 
exceptionally fine grids are used. The grids used for final results are as shown in Table 5.4. 


Table 5.4. Finite Element Model Discretization 


Panel Width, in. 

Total Number 
Of Elements 

Elements per 

Half-Wave 

X 

Y 

4.0 

620 

12 

8 

5.5 

528 

10/12 

8 

7.0 

456 

10.5 

10 
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A good starting point for grid evaluation is the bifurcation buckling analysis, since 
by examining convergence of the critical buckling load, one can get at least a preliminary 
estimate of the grid required. Shown in Table 5.5 are buckling loads obtained using several 
different grids in the analysis of the 4 inch- wide, [±30/90js panel. (The support elements 
around the edges of the model are omitted in the tabulation since they do not contribute 
to approximation of out-of-plane displacement.) Results shown in this table appear to be 
more sensitive to the number of elements used in the y direction than in the x direction; 
this is due to the fact that in the y-direction, the minimum number of elements used 
(4) is inadequate to approximate the out-of-plane behavior. For final results, the 12 by 
8 element model was chosen; however, based on limited experience with coarser grids, 
adequate results could be obtained with a somewhat lower level of discretization. 


Table 5.5. Grid Convergence, [±30/90]s 4-inch-wide Panel 


Elements/Half- Wave 
x-direction 

Elements/ Half-Wave 
y-direction 

Total 

Elements 

P 1 

* cr 

6 


192 

4142 

8 

4 

252 


12 


372 

4174 

6 



4247 

8 

8 


4272 

12 


620 

4290 

6 


448 

4278 

8 

12 

588 

4303 

12 


868 

4322 
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5.3.2 Solution Strategy 


The details of the nonlinear analysis procedure used herein are explained in terms of 
a typical end shortening curve for this type of problem, shown in Figure 5.7. The general 
idea is to allow STAGS to increment the load so that solutions are obtained at various 
points along the end shortening curve. The analysis is begun at a load below buckling, 
such as the point labeled one, and is continued far into postbuckling to a load several times 
the buckling load, such as point three. Since the analysis will be compared with test data, 
the final load is determined primarily by the panel’s exhibited failure load. 

The only difficulty is that, for the perfect case, the tangent stiffness matrix becomes 
singular and therefore unfavorable at buckling (point two in Figure 5.7). To overcome this 
problem, an imperfection is imposed on the panel in order to remove the singularity. The 
shape of the imperfection is obtained from a linear bifurcation buckling analysis, which 
generates the eigenvalue and mode shape and saves them on disk. Then, the user specifies 
that the nonlinear analysis will use the saved mode shape, multiplied by some constant. 
The result of this procedure is that, if the imperfection is not too large, the response 
is almost identical to that of the perfect case, except that the curve is differentiable at 
buckling, as shown by the solid line in Figure 5.7. The imperfection must not be too 
large, or the response will be grossly different from the perfect case. For the problems 
solved herein, the magnitude of the imperfection was about one percent of the total panel 
thickness, although a larger two percent imperfection was needed in the case of the seven 
inch wide, [±50/35] 5 panel. One might think that for panels with bending-extensional 
coupling, the out-of-plane behavior would be enough to trigger buckling without use of 
an imperfection, but such was not found to be the case. Without the imperfection, the 
analysis remains on the unstable primary path instead of branching to the stable secondary 
path, because the level of bending-extensional coupling present in the panels tested does 
not cause sufficient out-of-plane displacement to trigger buckling. 
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One more constraint is needed when the buckling analysis shows that the structure 
has two nearly coincident eigenvalues. In such a case, specifying the proper imperfection 
shape is not enough to cause the analysis to choose the correct secondary path; convergence 
difficulties arise and the analysis cannot get past buckling. In order to force the analysis 
to choose one of the two possible secondary paths, a single Lagrangian constraint can be 
applied to force the buckled shape to adhere to the symmetry of the desired mode shape. 
For example, if the desired mode is symmetric (meaning an odd number of half-waves), 
then the analyst specifies that two appropriate points on either side of the shorter centerline 
must have the same out-of-plane displacement; if the desired mode is antisymmetric, then 
the displacements at the two points must be negatives of one another. Of course, points 
must be chosen which represent the larger displacements in the panel, and the different 
types of symmetry exhibited by anisotropic panels (described in Chapter 2) must be kept 
in mind. 

The user must also choose the nonlinear solution strategy to be employed. STAGS 
provides essentially three solution techniques: the classical Newton-Raphson technique, 
a computationally less expensive modified Newton-Raphson technique, and an option in 
which the step size is controlled automatically using a path-length parameter as an inde- 
pendent parameter [24]. The path- length technique is used herein and was preferred over 
the two more classical Newton-Raphson techniques because they axe not as effective in the 
neighborhood of buckling, when the tangent stiffness matrix is ill-conditioned [25]. 

5.4 Comparison of Experimental and Finite Element Results 

5.4.1 Presentation and Discussion of Buckling Load 
And End Shortening Results 

There axe many aspects of the experiments which can be compared with computational 
results; perhaps the most informative of these in terms of global response is panel end 
shortening. This is because end shortening plots provide a simple measure of overall panel 
stiffness throughout the test, as well as a good estimate of the buckling load. There are 
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essentially three regions of interest in examining an end shortening curve. The first is 
the prebuckling range. It is important that agreement be obtained between model and 
experiment within this range, since agreement indicates that material stiffness coefficients 
and boundary conditions used in the analysis are correct. Secondly, one would like to 
predict the buckling load accurately. Finally, if the prebuckling and buckling responses 
axe accurate, then loading, boundary conditions, and material properties have probably 
been modeled correctly; therefore, any discrepancies in the postbuckling response cannot 
be due to improper modeling of these aspects of the problem unless they change as more 
load is applied. 

5. 4.1. a Buckling Load Results and Discussion 

A few comments are in order concerning the definition of the term “buckling load”. 
For a perfect panel with no bending-extensional coupling, the buckling load is very easy 
to determine, since buckling can be discerned as a sharp knee in the end shortening curve 
or the precise point at which back-to-back strain gages diverge. However, in the presence 
of initial imperfections or bending-extensional coupling, buckling is not well defined; in 
fact, some say that, to be precise, one should not even use the term “buckling” in such 
cases. This is because the transition from a prebuckled to a postbuckled configuration is no 
longer sudden; the phenomenon becomes gradual, and so the previously sharp bifurcation 
point in the end shortening and back-to-back strain gage curves becomes a gradual curve. 
However, in the interests of brevity, the term will be used herein, and will indicate that 
load at which a perfect system without bending-extensional coupling would buckle. 

The method used to estimate the buckling load is very simple and is depicted in 
Figure 5.7; on the end shortening curve, one draws two straight lines, one along the linear 
prebuckling portion of the curve, and one along the initial part of the postbuckling curve. 
The intersection of the two lines is then the buckling load, since an end shortening curve 
for the corresponding perfect system has been constructed. If desired, this estimated 
load can be confirmed by examining back-to-back strain data, as was shown in Figure 
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5.3. To evaluate this data, one examines the area of divergence of the two strains; the 
latter portions of the two curves are then extrapolated backward until they cross at some 
point within the region of divergence, as shown in Figure 5.3. This point is then taken as 
buckling. 

With this item of terminology defined, results may now be discussed. First, once the 
loading head leveling procedure described in Section 5.2.2 was used, excellent agreement 
was obtained for prebuckling stiffnesses, which were taken to be the slope of the end 
shortening curve in the linear prebuckling range. In all cases, the slope was predicted to 
within two percent or less. Secondly, buckling load results are shown in Table 5.6, where 
it is seen that average agreement between experiment and analysis was 14.4%. Both the 
experimental buckling loads, P®* p , and the STAGS buckling loads, P^ m , were obtained 
from the end shortening curves shown in Figure 5.8. The most likely explanation for the 
difference in experimental and analysis buckling loads is the presence of fairly large initial 
imperfections. Measurements indicated that imperfections of the order of 1.5 times the 
panel thickness were not uncommon. It should also be noted that the panel exhibiting the 
largest difference between analysis and experiment, the [±50/35]s, 4 inch-wide panel, also 
underwent a mode change, suggesting that the two phenomena might be related. 


Table 5.6. Buckling Loads, Pounds 


Panel Width, 
inches 

Panel Layup 

Number Of 
Half- Waves 

P c e r XP , lbs. 

P'V", lbs- 

Percent 

Difference 

4.0 

[±30/90] 5 

5 



14.7 

4. Of 

[±50/35] 5 

5— >6 



26.2 

5.5 


4 


mgm 

■ 

5.5 

[±50/35] s 

4 




7.0 

[±30/90] 5 

3 

mm 


12.0 

7.0$ 

[±50/35] s 

3 

wM 


7.9 


f Experienced mode change 
$ Required 2% imperfection 
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Linear bifurcation buckling analyses were also conducted for all the experiments dis- 
cussed; in all these analyses, the first two eigenvalues and mode shapes were calculated. 
In all but one case, the first mode determined by the analysis was also the mode exhib- 
ited by the panel. As shown in Table 5.7, the panel which was the exception underwent 
a mode change, with the second mode occurring first. This panel also had the smallest 
percent separation between its buckling loads, as well as the largest percent difference be- 
tween experimental and analysis buckling loads. As is obvious from the table, very close 
eigenvalues were frequent, but no trends suggest themselves regarding this phenomenon. 
It is also worth noting that all the panels exhibit the same trend that one would expect 
for isotropic or orthotropic panels: a higher aspect ratio means more half-waves in the 
buckling mode shape. 


Table 5.7. Linear Bifurcation Buckling Results, Pounds 


Width, 

Layup 

Mode 1 

Mode 2 

% Sep 4 

inches 

Load, lbs. 

Half-Waves 

Load, lbs. 

Half-Waves 

4.0 

[±30/90] 5 

4290 

5 

4311 

6 

■si 

4. Of 

[±50/35] 5 

5248 

6 

5265 

5 

m 

5.5 


3016 

4 

3028 

3 

0.4 

5.5 

BBS1 

3468 

4 

3559 

3 

2.6 

7.0 

[±30/90] 5 

2397 

3 

2458 

2 

2.5 

7.0 

[±50/35] 5 

2678 

3 

2839 

2 

5.7 


f Experiment exhibited mode change, with second mode occurring first; next eigenvalue 
is 5401 lbs. 


$ Indicates the percentage by which the buckling loads for the first two modes are 
separated. 


It is interesting to note that for these panels, comparison of Tables 5.6 and 5.7 indi- 
cates that the buckling loads obtained from the linear bifurcation buckling analyses are 
excellent estimates for the buckling loads obtained from the nonlinear analyses. The aver- 
age percent difference is only 1.5%. This result indicates that the prebuckling behavior is 
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not significantly nonlinear, since buckling loads are not appreciably affected by the use of 
nonlinear theory. 

5.4.1.b Postbuckling Results and Discussion 

The third and final portion of the end shortening curve depicts global postbuckling 
response. The end shortening curves shown in Figure 5.8 indicate varying degrees of 
agreement between experiment and analysis at two distinguishable stages of postbuckling. 
The first stage is a range of nearly linear behavior immediately after buckling. In this 
region, the slope of the end shortening curve suggests itself as a simple basis for comparison. 
Table 5.8 summarizes the slope values for the experiment and analysis curves, and indicates 
that agreement between analysis and experiment is good in this range, except in the case 
of the [±50/35] 5 , 4 inch- wide panel. Also, note that agreement is slightly better for the 
[±30/90] 5 panels, regardless of aspect ratio. 


Table 5.8. Slope of Initial Postbuckling Curve 


Layup 

Width, in. 

STAGS Slope 
SfxlO 4 

Experiment Slope 
£fxl 0 4 

Percent Error 


4.0 

5.89 

5.48 

7.5 

[±30/90] 5 

5.5 

7.47 

7.19 

3.9 


7.0 

9.43 

9.21 

2.4 


4.0 

2.77 

2.44 

13.5 

[±50/35] s 

5.5 

3.69 

3.51 

5.1 


7.0 

4.93 

4.58 

7.6 


The second discernible stage of postbuckling is marked by a gradual but noticeable 
degradation of the panel’s stiffness. Agreement between analysis and experiment also de- 
teriorates markedly in this range, with better correlation being obtained for the [±30/90] 5 
panels. By simply examining the end shortening plots, one can obtain a subjective esti- 
mate of the stiffness degradation of each panel; however, it is possible to obtain a numerical 
estimate which more clearly delineates the observed trends. Shown in Figure 5.9 is the 
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method used to obtain this estimate. In the figure, the early portion of the postbuckling 
curve is linearly extrapolated to a point having the same end shortening value as that 
of the experiment curve at a point far into postbuckling. Then, the difference between 
the extrapolated load and the actual load is taken as the stiffness degradation. The same 
procedure can be followed for the analysis data, and results can be compared; however, 
care must be taken to account for the differences in buckling load between experiment and 
analysis. To accomplish this, points on a set of experimental and analytical curves are 
chosen to correspond to the same load, nondimensionalized by the estimated buckling load 
P cr for the particular curve being examined. 

Results of these comparisons are shown in Table 5.9. Although comparisons are easily 
made between analysis and experiment, it is important to note that panel-to-panel com- 
parisons should not be made, since calculations were done at different values of P/P cr 
for each panel. Different values of P/P cr had to be chosen in order to obtain results 
which represented the stiffness degradation over the entire postbuckling range. Heeding 
this admonition, two trends are nevertheless apparent. First, the percent differences be- 
tween analysis and experiment are always significantly greater for the [±50/35]s panels. 
Secondly, for both layups, percent difference is largest for the 5.5 inch-wide panel, which 
was the only panel of the three to exhibit an antisymmetric buckling mode shape. 


Table 5.9. Postbuckling Stiffness Degradation 


Layup 

Width, in. 

P/Pcr 

STAGS, 

Ib/in 

Experiment, 

lb/in 

Percent 

Difference 


4.0 

2.0 

333 

400 

16.8 

[±30/90] 5 

5.5 

3.0 

800 

1267 

36.9 


7.0 

4.0 

2000 

2233 

10.4 


4.0 

1.4 

233 

333 


[±50/35] 5 

5.5 

2.0 

333 

867 



7.0 

2.5 

767 

1533 

50.0 
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The end shortening results presented here appear to be very similar to those presented 
in [16]. In [16], the effect of transverse shear on the various response quantities of rectangu- 
lar aluminum and [±45] graphite-epoxy panels is studied. It is important to note that the 
laminates and boundary conditions used in [16] are different from those described herein; 
thus definite conclusions are difficult to draw based on comparison with the present work. 
Nonetheless, some observations can still be made. 

The end shortening curves presented in [16] show that, for the cases studied, transverse 
shear theory predicts a slightly less stiff postbuckling response than does classical theory. 
These results are similar to those shown in Figure 5.8, in that the experimental results are 
less stiff than the finite element results obtained without shear deformation. Thus, it is 
tempting to attribute the lack of test /analysis correlation dining postbuckling to the lack of 
a transverse shear flexible finite element formulation. However, it is important to note that, 
in [16], the differences between the classical and transverse shear results begin immediately 
at buckling; there is no postbuckling region in which the stiffness remains nearly constant 
and the slopes of the two curves agree, as is the case with the results presented in Table 5.8. 
Furthermore, the results presented in Table 5.9 differ greatly from those observed in the 
graphs of [16]. There is no significant disparity in the stiffness degradation predicted by 
the shear deformation and classical theories in [16], while the disparity is quite large for the 
panels examined in the present work. For these two reasons, it is difficult to attribute all 
the differences between experiment and finite element analysis to a lack of transverse shear 
flexibility in the finite element model; however, the role of transverse shear in the response 
of anisotropic laminates with bending-extensional coupling should not be discounted, and 
merits further study. 

A second explanation for the degradation of postbuckling stiffness in all the experi- 
ments is that a material failure of the type described in [26] has occurred. In [26], the point 
is made that because of the form of the stress transformation relations for an all-±45 lami- 
nate, large shear stresses are encountered in the principal material directions and therefore 
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at the matrix/fiber interface of each lamina. These shear stresses tend to break the bond 
between the fiber and matrix, resulting in an inplane matrix shearing failure. This failure 
mechanism seems a likely cause for the more severe postbuckling stiffness degradation ob- 
served in the [±50/35]s laminate, as well as for the larger disparity between predicted and 
observed stiffness degradation, since this laminate is very similar to an all- ±45 laminate. 
Also, some of the C-scans shown in Appendix C support the existence of matrix shearing, 
since there does appear to be some damage which runs in the fiber direction. 

5.4.2 Presentation and Discussion of Out-of-Plane Displacement Results 

Out-of-plane (or transverse) displacement data was collected at seven to eight differ- 
ent points on each of the panels tested. Data obtained at the point of maximum trans- 
verse displacement is of interest, since it allows determination of the relative magnitude 
of out-of-plane behavior at the different stages of the experiment. Also, as mentioned 
previously, comparison of transverse displacement at reflectionally symmetric points in- 
dicates the amount of twist occurring throughout the experiment. Lastly, satisfaction of 
the anisotropic symmetry conditions described in Chapter 2 can be checked by comparing 
results at inversionally symmetric points. 

Shown in Figure 5.10 are plots of maximum transverse displacement for all six panels. 
From these plots, three facts can be noted. First, the maximum transverse displacements 
occurring in these experiments range between one and four times the panel thickness 
of approximately 0.08 inches. Second, the transverse displacement present at buckling 
amounts to less than 40% of the panel thickness, decreases with increasing aspect ratio, 
and is generally larger for the [±30/90]s panels, as summarized in Table 5.10. Finally, 
the STAGS curves and the experimental curves have almost exactly the same curvature in 
postbuckling, indicating that if one adjusts for the error in buckling load, the transverse 
displacements are modeled accurately. 

The plots of Figure 5.11 indicate the amount of twisting occurring in each of the panels. 
The two curves appearing in each plot are the absolute values of transverse displacements 
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Table 5.10. Transverse Displacement at Buckling 


Layup 

Width, in. 

Displacement, in. 

Displacement 

Thickness 


4.0 

0.020 

26 

[±30/90] 5 

5.5 

0.025 

32 


7.0 

0.029 

37 


4.0 

0.006 

8 

[±50/35] 5 

5.5 

0.017 

22 


7.0 

0.016 

21 


at two reflectionally symmetric points on either side of the short centerline, as shown in the 
figure legend. If the two lines are coincident, then no twisting has occurred; the distance 
between the two lines is then the amount of twist. In the case of the [±50/35]s, 4 inch- wide 
panel, a mode change prevents meaningful comparison of the two displacements; therefore, 
discussion is presented here for the remaining five plots. 

Overall, there is no pattern to the presence or absence of twisting behavior. Notice 
that for the [±30/90] 5 panels, little or no twist is observed when the buckling mode is 
symmetric (i.e., for the 4 and 7 inch-wide panels), but a significant amount of twist occurs 
in the 5.5 inch panel, which experienced an antisymmetric buckling mode. Prom these 
results, one might be tempted to assert that twisting occurs only when the buckling mode 
is antisymmetric. However, for the remaining two [±50/35]s panels, the amount of twist 
is non-negligible regardless of mode shape, thus destroying the correlation between twist 
and mode shape symmetry. Shown in Figure 5.12 are corresponding twist plots generated 
from the finite element analysis results. These plots are consistent with Figure 5.11. 

The plots shown in Figure 5.13 are similar to the twist plots, except that they are 
displacement plots at inversionally symmetric points, and so indicate the degree of symme- 
try satisfaction. As was explained in Chapter 2, no precise type of symmetry is expected 
when an antisymmetric buckling mode is exhibited, and this is confirmed in the plots for 
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both 5.5 inch-wide panels. For the [±30/90]s, 4 inch-wide panel, the symmetry violation 
is of the same magnitude as the twist, i.e., very small, while for the 7 inch-wide panels, 
the symmetry conditions are fairly well satisfied early in postbuckling, but become badly 
violated as the experiment proceeds. Shown in Figure 5.14 are similar plots made from 
the STAGS results; these plots confirm all the symmetry conditions stated in Chapter two, 
including the fact that there is no symmetry when the buckling mode is antisymmetric, as 
in the case of the 5.5 inch-wide panels. 

To summarize the twist and symmetry results, the following table is presented: 


Table 5.11. Summary of Twist and Symmetry Results 


Layup 

Width, 

Symmetry Of 

Twist Magnitude 

Symmetries Satisfied? 

inches 

Mode Shape 

Experiment 

STAGS 

Experiment 

STAGS 

[±30/90] 5 

4 

7 

Symmetric 

Negligible 

Zero 

Yes 

No 

Yes 


5.5 

Antisymmetric 

Large 

Large 

None Expected or Observed 


4 

— 

— 

— 

— 

— 

[±50/35] 5 

7 

Symmetric 

Large 

Large 

No 

Yes 


5.5 

Antisymmetric 

Large 

Large 

None Expected or Observed 


5.4.3 Presentation and Discussion of Strain Results 

Back-to-back strain gage rosettes and axial strain gages were placed at between nine 
and fifteen locations on the tested panels. Attempts were made in all cases to record 
maximum or near-maximum strains, and in two cases, strains were recorded near failure 
locations. 

Shown in Figure 5.15 are plots of strain at those points on each panel which underwent 
maximum and near maximum strains. In the case of the 5.5 inch-wide panels, the maximum 
strain, which occurred in the center of the second of four half-waves, was not recorded; 
instead, strain in the center of the first half-wave was measured. For the other four panels, 
center strains were the maximums, since an odd number of half-waves were observed in 
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those cases. Also shown in Figure 5.15 are the estimated buckling loads for each of the 
panels, which are confirmed by the plotted strains. 

The most notable characteristic of these six strain plots is the unusual shape of most of 
the curves. Look first at plot (b), the center strains for the [±50/35]s, 4 inch-wide panel. 
The discontinuities late in postbuckling are easily explained; they indicate an observed 
mode change from a five to a six half-wave buckle pattern, and so axe not unexpected. 

Plots (c) and (d), the near-maximum strains for both 5.5 inch-wide panels, are very 
similar to other published back-to-back strain curves [15], in that the two strain curves 
grow smoothly apart and do not curve back toward each other at any point. Since the 
distance between the curves is essentially twice the local curvature in the axial direction, 
such behavior simply indicates that the initial buckled shape is becoming more pronounced; 
i.e., the curvature is steadily increasing. Now compare these plots to plot (f). There, some 
curvature is developed immediately after buckling, and then about one-third to one-half 
of the way through postbuckling, the curvature decreases. Such behavior is often taken to 
indicate an overall mode change; however, this seems peculiar since the plot of maximum 
transverse displacement indicates that the maximum out-of-plane deformation increases 
smoothly throughout postbuckling. Also, no mode change was observed in the moire 
fringe photographs. At most other points of non-negligible strain on the same panel, the 
back-to-back strain curves are either very similar to those shown in Figure 5.15, or the 
two strain curves actually cross over each other, indicating that the curvature has changed 
sign. 

The most obvious explanation which accounts for both the displacement and strain 
behavior is that, in a gross sense, the initial buckled shape is retained throughout post- 
buckling, while some of the details of the deformation pattern change. Figure 5.16, a 
plot of transverse displacement down the long centerline of the [±50/35]s, 7 inch-wide 
panel, displays just such behavior. (This plot was generated using STAGS results, and the 
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shown strains are actual test data.) Initially, the deformed shape appears nearly sinusoidal. 
However, at P = 5733 pounds, the following changes have taken place: 

• The center half-wave has flattened, causing a decrease in curvature as shown in plot 
(b) of the figure, while the outer two half-waves have become more pointed, explaining 
the increase in curvature displayed in plot (a). 

• The points of zero deflection (commonly referred to as the nodal lines) have shifted 
substantially. This explains the curvature sign change seen in plot (c); this set of 
gages was placed near enough to the nodal line for the fine to have shifted from a 
position above the strain gages to a point below them. No sign change occurred for 
gages placed far away from a nodal line. 

• The inflection points, denoted by the different symbols, have shifted, indicating that 
the entire deflection pattern has shifted off the load axis. For brevity, this phenomenon 
will henceforth be referred to as pattern eccentricity. (The inflection points were 
determined using cubic spline interpolation of the STAGS displacement data.) 

Finally, when the load has increased to 8000 lb, curvature at the center of the panel has 
increased slightly, thus explaining the “hourglass” shape of plot (b). 

Similar plots made for the other two [±50/35 ] 5 panels display all the same trends 
as Figure 5.16, except that for the 4 inch-wide panel, pattern eccentricity is negligible. 
(See Figure 5.18 for comparison.) On the other hand, plots produced from the [±30/90]s 
panel analysis data behave somewhat differently, as shown in Figure 5.17. Comparison of 
this figure and Figure 5.16 indicate that the flattening effect is much more severe in the 
[±50/35]5 panel; in the [±30/90]s panel, only enough flattening occurs to cause a barely 
noticeable curvature change in strain plot (b). A shift in the location of the nodal lines 
is still observed for this panel, causing a curvature sign change in plot (c); also, the outer 
half-waves become progressively more pointed with load, again causing the steady increase 
in curvature depicted in plot (a). Pattern eccentricity is not seen in this figure. 
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Further examination of all the plots shown in Figure 5.18 reveals that one general- 
ization is possible: the flattening effect is consistently more exaggerated in the [±50/35] 5 
panels. For the [±30/90 ] 5 panels, curvature at points of maximum displacement either 
remains nearly constant after a certain load, or it increases. The other two characteris- 
tics, shift of nodal lines and pattern eccentricity, occur with no particular pattern; that 
is, no real generalizations regarding these phenomena can be made. Furthermore, both 
phenomena can be observed in orthotropic panels. 

It is apparent from this discussion that decreases in curvature partway through post- 
buckling do not necessarily indicate a mode change, at least not in the classical sense of a 
shift from one deformation pattern to another, dramatically different pattern because of a 
newly encountered instability. Rather, such curvature changes can point to subtle changes 
in the deformation pattern which can best be described as changes in the overall character 
of the buckled mode shape. 

It was suggested by Jensen [14] that such curvature changes represent a localized 
dimpling, indicating the onset of a classical mode change which actually did occur for his 
panels. Such a phenomenon is possible and could account for the localized presence of 
unusual strain patterns; however, the presence of such patterns at many points on the 
plate, as well as the consistent analysis and test displacement data presented here suggest 
that dimpling is not causing the decreases in curvature; rather, gradual flattening of central 
portions of the pattern is at fault. 

5.5 Failure Initiation and Development 
5.5.1 Qualitative Results 

From the combination of the moire fringe photographs (shown in Appendix B), the 
C-scan results (shown in Appendix C), and the failed panels themselves, a good qualitative 
picture of the failure characteristics of each panel is obtained. Shown in Figure 5.19 axe 
sketches, drawn to scale, of all six panels. Three things are indicated on each sketch: area 
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of failure initiation as detected by C-scans, area of final failure, and final position of the 
nodal lines as indicated by the moire fringe photographs. No measurements are shown, 
since the photographs do not allow precise measurement. Determining the locations of the 
nodal lines is particularly difficult, since it is hard to ascertain the exact locations of lines 
of zero displacement by measuring the photographs. 

In any case, a few observations are possible. First, note that failure always initiates 
either on or very close to a nodal line. This fact is consistent with [15]. There, the ex- 
planation is advanced that “the higher strains near the specimen edges couple with the 
out-of-plane deflection gradients at the nodal line to induce sufficient transverse shearing 
loads to fail the specimen in shear before the large bending strains at the points of max- 
imum displacement become critical”. The facts embodied in this statement can for the 
most part be confirmed by information presented in various parts of this chapter, as will 
be explained in the following paragraph. 

First, it is easy to confirm that strains near the edges of the panel are higher than 
those along the centerline. Shown in Figure 5.20 are plots of strain data collected at five 
closely spaced points near the short edges of both 5.5 inch-wide panels. This data allows 
comparison of strains over the width of both panels. The magnitudes of the middle surface 
strains, indicated by open circles, are always higher near the edges of the panels. (Note, 
however, that this is not always true of the surface strains; near the bottom of the [±30/90] 5 
panel, one of the center surface strains is actually slightly larger than the two edge strains.) 
Secondly, strain data collected near the nodal lines of the 7 inch-wide panels confirms that 
the strains are much higher along the specimen nodal lines, as shown in Figures 5.16 and 
5.17. (Whether strain is higher at the edge of the nodal line than at its center cannot be 
determined here, since no data was taken at the center of a nodal line.) Third, it is obvious 
from these same two figures that the gradient of transverse displacement is also high at the 
nodal line. Finally, that the panels fail in shear is suggested by the photograph shown in 
Figure 5.21. To obtain this photograph, the [±50/35]s, 7 inch-wide panel was sectioned so 
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that the area where the C-scan indicated failure initiation could be examined. A diagonal 
crack through the thickness is evident, with delaminations propagating in both directions 
away from the crack. 

A second observation involves the relationship between initial and final failure loca- 
tions. In four of the six panels, the initial delamination and the final failure were nowhere 
near each other. This is especially true for the 7 inch-wide panels, where the C-scans 
indicated failure initiation at the top nodal line, while the final failure occurred at the 
bottom nodal line. Evidently, some redistribution of load occurs after the initial damage is 
sustained, causing the final failure to frequently occur in a totally different location from 
the initial damage. 

5.5.2 Quantitative Results 

Table 5.12 presents a summary of pertinent quantitative data collected at the failure 
load of each panel. From this data, the following facts are apparent: 

• Normalized loads increase steadily with decreasing aspect ratio, as expected. However, 
absolute failure loads are nearly the same for the 4 and 5.5 inch-wide panels of each 
layup, while the load withstood by the 7 inch-wide panels is somewhat higher for both 
layups. 

• Loads carried by the [±30/90] 5 panels are much higher, both in absolute and relative 
terms, than those carried by the [±50/35]s panels. 

• Normalized end shortening increases with decreasing aspect ratio, while absolute end 
shortening does just the opposite. 

5.6 Summary of Observations 

The following is a short summary of observations made in this chapter. Refer to 
Chapter 6 for discussion of conclusions, as well as for suggestions for future research. For 
the reader’s convenience, notation is made below of the sections, tables, and figures which 
best support each item. 
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Table 5.12. Load and End Shortening at Failure 


Layup 

Width, in. 

Pf, ^s. 

Pf/Pcr 

- - ... 

6/, in. 

bf/^cr 


4.0 

8500 

2.27 

.129 

3.9 

[±30/90] 5 

5.5 

8570 

3.12 

.119 

6.6 


7.0 

9480 

4.31 

.112 

9.7 


4.0 

6790 

1.62 

.225 

2.7 

[±50/35]5 

5.5 

6740 

2.17 

.199 

4.7 


7.0 

7220 

2.93 

.195 

7.4 


Observations based on experiment behavior: 

• The experimental setup is critical for laminates possessing anisotropy and/or bending- 
extensional coupling, in that mirror symmetry conditions cannot be used to ensure 
that a uniform axial load is applied. (Section 5.2.1.) 

• Buckling loads tend to be slightly higher for the [±50/35)5 panels, with the differences 
between results for the two layups being larger for the higher aspect ratio panels. 
(Table 5.6.) 

• Transverse displacements and bending strains at buckling are not appreciable, in- 
dicating that the effect of bending-extensional coupling on prebuckling behavior is 
negligible for the laminates studied. (Figure 5.15, Table 5.10, Tables 5.6 and 5.7, 
Section 5.4.1.) 

• According to end shortening results, there are two discernible regions in postbuckling. 
During the initial portion of postbuckling, the initial postbuckling stiffness of the panel 
is maintained. Then, about halfway through the postbuckling range, the stiffness of 
the panel begins to drop substantially. (Figure 5.8) 

• Strain results indicate subtle changes in the buckle mode shape throughout postbuck- 
ling. The most noticeable changes are dramatic decreases in curvature in the central 
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portions of the [±50/35]s buckle patterns, as well as pronounced shifting of nodal lines 
for some panels. (Section 5.4.3.) 

• Initiation of failure always occurs at or very close to a nodal line of the buckled shape. 
This lends credence to the observations advanced in [15]. (Figure 5.19.) 

• Locations of failure initiation and final failure are not necessarily the same, probably 
due to redistribution of load after the initial damage. (Figure 5.19.) 

• Postbuckling load carrying capability is much greater for the [±30/90 ] 5 panels than 
for the [±50/35]5 panels. (Table 5.12.) 

• Maximum observed transverse displacement ranges between one and four times the 
panel thickness and is inversely proportional to panel aspect ratio. (Figure 5.10.) 

• Satisfaction of displacement symmetry conditions described in Chapter 2 is not con- 
sistent; the conditions are violated in both of the 7 inch-wide panels. (Section 5.4.2.) 

• Twisting behavior is observed for both layups, but does not appear in all six panels. No 
pattern of absence or presence of twisting is discernible, since no twisting is observed, 
either analytically or experimentally, in the two higher aspect ratio, [±30/90]s panels. 
(Section 5.4.2.) 

Observations based on finite element modeling efforts: 

• Modeling of global response in the prebuckling range presents no problem, since the 
test and analysis end shortening results differed less than two percent within the 
prebuckling range. (Section 5.4.1.) 

• Buckling loads estimated from a nonlinear finite element analysis are an average of 
14.4% different from buckling loads estimated from test results. Predictions are better 
for lower aspect ratio panels, and are slightly better for [±30/90]5 panels (13.1% 
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average difference) than for [±50/35]s panels (15.7% average difference). The most 
likely cause of the overall disparity is the presence of imperfections of the order of 1.5 
times the panel thickness. (Table 5.6.) 

• For the panels examined herein, linear bifurcation buckling analysis is as accurate as 
a nonlinear analysis in the prediction of buckling loads. This result again emphasizes 
that the prebuckling behavior of the panels studied was not appreciably nonlinear. 
Bifurcation analyses were also very helpful in determining expected mode shapes, as 
well as in pointing out the likelihood of a mode change during postbuckling. (Tables 
5.6 and 5.7.) 

• Linear bifurcation buckling analysis predicts nearly coincident eigenvalues, i.e., eigen- 
values separated by less than .5%, for three of the six panels. Only in the case where 
separation of the eigenvalues was the smallest (.3%) was an obvious mode change 
observed. (Table 5.7.) 

• In the initial portion of postbuckling, agreement between test and analysis end short- 
ening results is good (2% to 8% difference in slope) except in the case of the [±50/35] 5 , 
4 inch- wide panel, where the slopes differ by 13.5%. In all cases, slightly better agree- 
ment is obtained for the three [±30/90]s panels. (Table 5.8.) 

• In the later portion of postbuckling, the stiffness degradation is always larger than that 
predicted by finite element analysis. Predictions of the drop in stiffness (as measured 
by percent difference between analysis and test) are much better for the [±30/90] 5 
panels by factors of two to five, and are significantly worse for the 5.5 inch-wide panels 
of both layups. Possible causes for this disparity between analysis and test include the 
lack of transverse shear flexibility in STAGS and the possibility that matrix shearing 
occurs long before final failure. (Table 5.9, Section 5.4.1.) 
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• Finite element results indicate subtle changes in the buckling mode shape occurring 
throughout postbuckling. In all the [±50/35]s panels, the mode shape tends to flat- 
ten, causing decreases in curvature at points of maximum transverse displacement. 
This observation is confirmed by strain results, as described above. In some panels, 
pronounced shifting of nodal lines and inflection points also occurs, but no pattern is 
evident in this behavior, and such phenomena can be observed in orthotropic panels. 
(Section 5.4.3.) 

• In all three cases, the [±50/35]s panels have proven more difficult to model in post- 
buckling than the [±30/90]s panels. In reviewing the results for the two sets of panels, 
the fact stands out that the largest differences between analysis and experimental re- 
sults were usually seen for the [±50/35]s panels. In particular, difficulties arose in 
predicting the slope of the initial part of the postbuckling end shortening curve, and 
the stiffness degradation during the latter part of postbuckling. The only substantial 
difference between the two panels is that the [rt30/90]s panel does not possess any 
extension-shear coupling; however, the evidence is not sufficient at this time to at- 
tribute the modeling difficulties to this distinguishing feature. (Throughout chapter.) 

• The one panel which exhibited an observable mode change ([±50/35] 5 , 4 inch-wide) 
has proven very difficult to model. Apart from the current lack of a simple way to 
model the actual mode change, other seemingly unrelated aspects of the behavior were 
more difficult to model. For example, the percent difference in buckling load between 
analysis and experiment was 12.5% greater for this panel than for any of the others. 
Similarly, initial slope of the postbuckling curve was predicted 5.9% worse for this 
panel than for any other. (Throughout chapter.) 

• Transverse displacements are well-predicted if adjustments are made to account for 
the error in predicting the buckling load. (Figure 5.10.) 
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• Although symmetry satisfaction is inconsistent in the experimental results, STAGS 
results behave precisely as predicted by the symmetry conditions of Chapter 2. (Table 
5.11.) 
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Chapter 6 

Observations and Recommendations 


6.1 Observations 

Following is a summary of the observations derived from the present work. For the 
reader’s convenience, notation is made of the section, tables, and figures which best support 
each statement. Also, each point is preceded by a short statement which categorizes it 
within the subsection. 

6.1.1 Observations Related to the Modified Ravleigh-Ritz Technique 
6. 1.1. a Observations related to linear results 

• Effect of boundary conditions on results: Accurate solution of problems involving 

the two axially loaded panels required the use of more functions per displacement 
degree of freedom than did similar problems involving transverse loading. This was 
due to the more complex boundary conditions used in the axial loading problems; in 
those problems, the in-plane displacement u 2 was constrained along the loaded edges, 
causing a complicated stress state to be present in the corners of the panels. (Sections 
4.2.3 and 4.2.4) 

• Effect of anisotropy and bending-extensional coupling on results: Accurate solution 

of problems involving the anisotropic panel with bending-extensional coupling required 
the use of many more functions per displacement degree of freedom than did identical 
problems involving an orthotropic panel. (Section 4.2) 

• Effectiveness of exterior penalty approach for enforcing boundary conditions: 

The external penalty function approach was effective in assuring satisfaction of bound- 
ary conditions for all problems, although in the case of the axially loaded panels, 
satisfaction of the constraint 112 = 0 constraint along the loaded edge was somewhat 
sensitive to variations in the penalty parameter. Overall, a broad range of four orders 
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of magnitude in the penalty parameter gives acceptable results for a linear problem; 
however, a narrower range of two orders of magnitude is recommended for nonlin- 
ear analysis due to observed convergence problems attributable to too high a penalty 
parameter. (Section 4.3) 

• Performance of basis vectors: 

o For the problems involving axial loading of anisotropic panels, the chosen basis vectors 
performed well. (Section 4.4.1) 

o For the problems solved herein, the Grammian provides a good estimate of the number 
of basis vectors required to solve a given problem, thus alleviating the need to examine 
detailed results in order to select an appropriate reduced system size. (Section 4.4.2) 

6.1.1.b Observations related to nonlinear results 

• Effect of anisotropy and bending-extensional coupling on results: Using the modi- 

fied Rayleigh-Ritz technique, postbuckling response of the orthotropic panels studied 
can be accurately predicted; however, about twice as many functions are required in 
order to predict the postbuckling response of the nonorthotropic panels studied. A 
reasonable estimate of the postbuckling behavior of the nonorthotropic panels was 
obtained using the largest possible number of approximation functions (as defined by 
computer hardware limitations.) (Section 4.5) 

• Effect of first order transverse shear deformation effects: Slight differences in over- 

all postbuckling stiffness were seen which could definitely be attributed to transverse 
shear deformation, as in the case of the [04/904] s panel; however, the addition of trans- 
verse shear deformation did not result in better prediction of postbuckling stiffness 
degradation in either of the other two cases studied. (Section 4.5) 
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6.1.2 Observations Based on Experiments 

• Experimental setup: The experimental setup is critical for laminates possessing 

anisotropy and/or bending-extensional coupling, in that mirror symmetry conditions 
cannot be used to ensure that a uniform axial load is applied. (Section 5.2.1.) 

• Prebuckling range: Transverse displacements and bending strains at buckling are not 
appreciable, indicating that the effect of bending-extensional coupling on prebuckling 
behavior is negligible for the laminates studied. (Figure 5.15, Table 5.10, Tables 5.6 
and 5.7, Section 5.4.1.) 

• Postbuckling: 

o According to end shortening results, there are two discernible regions in postbuckling. 
Dining the initial portion of postbuckling, the initial postbuckling stiffness of the panel 
is maintained. Then, about halfway through the postbuckling range, the stiffness of 
the panel gradually begins to drop. (Figure 5.8) 

o Strain results indicate subtle changes in the buckle mode shape throughout postbuck- 
ling. The most noticeable changes are dramatic decreases in curvature in the central 
portions of the [±50/35] 5 buckle patterns, as well as pronounced shifting of nodal lines 
for some panels. (Section 5.4.3.) 

o Experimental satisfaction of displacement symmetry conditions described in Chapter 
2 is not consistent. (Section 5.4.2.) 

o Twisting behavior is observed for both layups, but does not appear in all six panels. 
No pattern of absence or presence of twisting is discernible. (Section 5.4.2.) 

• Failure: 

o Initiation of failure always occurs at or very close to a nodal line of the buckled shape. 
This lends credence to the observations advanced in [15]. (Figure 5.19.) 

o Locations of failure initiation and final failure are not necessarily the same, probably 
due to redistribution of load after the initial damage. (Figure 5.19.) 
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6.1.3 Observations Based on Finite Element Modeling Efforts 


• General observations: 

o In all three cases, the experimental postbuckling behavior of the [±50/35] 5 panels was 
more difficult to predict analytically than that of the [±30/90] 5 panels. In reviewing 
the results for the two sets of panels, the fact stands out that the largest differences 
between analysis and experimental results were usually seen for the [±50/35] 5 panels. 
In particular, difficulties arose in predicting the slope of the initial part of the post- 
buckling end shortening curve, and the stiffness degradation during the latter part 
of postbuckling. The only substantial difference between the two panels is that the 
[±30/90] 5 panel does not possess any extension- shear coupling; however, the evidence 
is not sufficient at this time to attribute the difficulties to this distinguishing feature. 
(Throughout chapter 5.) 

o The experimental behavior of the one panel which exhibited an observable mode 
change ([±50/35] 5, 4 inch- wide panel) has proven very difficult to predict analytically. 
Apart from the current lack of a simple way to model the actual mode change, other 
seemingly unrelated aspects of the behavior were more difficult to model. For example, 
the percent difference in buckling load between analysis and experiment was 12.5% 
greater for this panel than for any of the others. Similarly, prediction of initial slope of 
the postbuckling curve was 5.9% worse for this panel than for any other. (Throughout 
Chapter 5.) 

• Prebuckling range: Modeling of global response in the prebuckling range presents 

no problem, since the test and analysis end shortening results differed less than two 
percent within the prebuckling range for all panels. (Section 5.4.1.) 

• Buckling: 

o Buckling loads estimated from a nonlinear finite element analysis are an average of 
14.4% different from buckling loads estimated from test results. A likely cause of the 
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overall disparity is the presence of imperfections of the order of 1.5 times the panel 
thickness. (Table 5.G.) 

o For the panels examined herein, linear bifurcation buckling analysis as obtained from 
STAGSC-1 is as accurate as a nonlinear analysis in the prediction of buckling loads. 
Bifurcation analyses were also very helpful in determining expected mode shapes, as 
well as in pointing out the likelihood of a mode change during postbuckling. (Tables 
5.6 and 5.7.) 

• Postbuckling: 

o In the initial portion of postbuckling, agreement between test and analysis end short- 
ening results is good for five of the six panels. (Table 5.8.) 

o In the later portion of postbuckling, the stiffness degradation is always larger than that 
predicted by finite element analysis. Predictions of the drop in stiffness (as measured 
by percent difference between analysis and test) are much better for the [±30/90] 5 
panels by factors of two to five, and are significantly worse for the 5.5 inch-wide panels 
of both layups. Possible causes for this disparity between analysis and test include the 
lack of transverse shear flexibility in STAGS and the possibility that matrix shearing 
occurs long before final failure. (Table 5.9, Section 5.4.1.) 

o Finite element results indicate subtle changes in the buckling mode shape occurring 
throughout postbuckling. In all the [±50/35 ] 5 panels, the mode shape tends to flatten, 
causing decreases in curvature at points of maximum transverse displacement. This 
observation is confirmed by strain results, as described in the previous section. In 
some panels, pronounced shifting of nodal lines and mode shape inflection points 
also occurs, but no pattern is evident in this behavior, and such phenomena can be 
observed in orthotropic panels. (Section 5.4.3.) 
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o Satisfaction of displacement symmetry conditions described in Chapter 2 is consistent 
for the finite element results; the conditions are violated only for the two panels which 
exhibited antisymmetric buckle mode shapes. (Section 5.4.2.) 

6.2 Recommendations 

6.2.1 Recommendations Related to the Modified Ravleigh-Ritz Technique 

• Work should be done to obtain a more suitable set of approximation functions for 
analysis of anisotropic panels, continuing to bear in mind the symmetry relations 
for such panels which could be utilized to simplify the selection process. Using the 
existing Rayleigh-Ritz code, new sets of functions can be easily investigated. 

• A more efficient nonlinear solution technique such as an arc-length projection tech- 
nique which would work in conjunction with the existing Newton-Raphson technique 
should be implemented in order to minimize program execution times. 

• The reduced basis techniques discussed in Section 3.4. l.b should be implemented using 
the existing code. The code was written with an eye toward doing this, including 
organized generation of all the necessary partitioned arrays utilized in Section 3. 4. l.b. 

6.2.2 Recommendations Related to Experimental Work 

• Additional studies should be conducted to determine the cause of the gradual degra- 
dation of stiffness late in the postbuckling range. Panels with mechanical couplings 
similar to those used by Jensen in [14] should be tested in the same way as those exam- 
ined herein in order to isolate the cause of the phenomenon. In particular, laminates 
which are likely to display matrix shearing failures (e.g., ±45-dominated laminates) 
should be carefully isolated, since matrix shearing is one possible cause of the observed 
stiffness degradation. 

• Along a similar vein, additional attention should be focused on determining precisely 
why the [±50/35]s panels were generally more difficult to model than the [±30/90]s 
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panels. Again, laminates which axe likely to exhibit matrix shearing should be treated 
with care. 

• More tests involving panels which exhibit mode changes should be performed, since 
this behavior appears to introduce seemingly unrelated modeling difficulties into the 
problem. Such tests should rely on the panel geometry to induce mode changes, rather 
than bringing in the unnecessary complications of anisotropy and bending-extensional 
coupling. Apparently, simple bifurcation buckling analyses using STAGS can provide 
some idea of which panels are likely to exhibit mode changes, since the panel which 
displayed a mode change also had eigenvalues spaced only 0.3% apart, the smallest 
spacing of any of the six panels. 

• Additional tests involving anisotropic panels should be performed in which severed 
DCDT’s (perhaps 10) are placed down the centerline of the panels in order to further 
study the gradual flattening of the postbuckled shape with increases in load. Or- 
thotropic panels should also be tested in this way as control specimens, to confirm 
that the phenomenon is due to anisotropy. 
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Appendix A 

Array Symmetrization Using Only Quadratic Energy Terms 


This appendix will outline briefly the steps which must be taken in order to “sym- 
metrize” the stiffness array Kij arising from the quadratic strain energy term 


u 2 = \k ij x i x j 


(A.l) 


Note that Kij is already symmetric; nevertheless, symmetrization can still be used, if only 
for illustration purposes. 

In order to take the variation of U2, differentiation must be carried out with respect 
to X a , where a is a dummy index having the same range as I and J. Therefore, 

SU 2 = ( A - 2 ) 

aX a 

The product rule of differentiation is then applied, as follows: 

w ’-Mi: Xj+x @“* (A3) 


Since dXj/dX a can be written as 6j a , and similarly, dXj/dX a = 6j a , where 

*> = { 0 


l =J 

* 7^ J 


(A. 4 ) 


the variation of U2 may be rewritten as 

6U2 — -KuiSiaXj + 'X.i6j a )6X a 

= \ (K a jXj + K Ia X f ) 6X a (A. 5 ) 

= ±(K I j + K ji )Xj6X i 

Notice that the last equation involves the sum of a matrix and its transpose, divided 
by two, which is the definition of the symmetric part of the matrix Kij and will be 
denoted by the unbarred symbol Kjj. This same kind of “symmetrization” also occurs as 
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a natural consequence of the variational procedure for both of the nonlinear arrays Fijk 
and Gukl- For Fi jk, the symmetrized array Fjjk consists of a factor of £ and the sum 
of six permutations of Fijk > while for Gijkl , a factor of ^ and twenty-four permutations 
of Gijkl are required. Note that in the case of the three- and four- dimensional arrays 
Fijk and Gijkl-, symmetry of these arrays means that interchanging any two indices does 
not change the array. 

Finally, the incorporation of the symmetrized array K into the expression for 8U2 
gives 

8U 2 = KijXjSXi (A. 6 ) 

The above result, without £Xj, can be thought of as contributing to the right hand side of 
the Newton-Raphson equations. Obtaining the contribution to the left hand side is trivial; 
it is simply K. 

Again note that this symmetrization procedure is not particularly important for the 
quadratic energy term since the matrix K ij was symmetric to begin with; however, for the 
cubic and quartic energy terms which give the Fijk and Gijkl arrays, symmetrization 
is of paramount importance, since those two arrays are not symmetric before the variation 
is taken. In implementing a computerized nonlinear solution technique involving these 
arrays, a significant storage savings can be realized if the symmetry of the arrays is noted, 
since all three can become quite large. 
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Appendix B 

Experimental Moire Fringe Photographs 

This appendix contains representative Moire fringe photographs for each of the six 
tests discussed in Chapter 5. Each photograph includes the load meter, which displays the 
applied load in thousands of pounds. 

For each panel, four photographs are presented and are labeled as follows: 

1. Before Test: photograph of the test setup prior to introduction of load. This 
photograph gives a qualitative indication of the initial imperfection present in 
the panel. 

2. Buckling: photograph taken at a load as close as possible to the experimental 
buckling load as presented in Chapter 5. 

3. Late Postbuckling: photograph taken at a load far into the postbuckling range. 

4. Failure: photograph taken after the panel has failed. The load shown in this 
photograph will not correspond to the failure load reported in Chapter 5, since it 
was impossible to capture the panel at the precise moment of failure, before the 
panel’s load carrying capability dropped dramatically. 

One additional set of two photographs is included for the [±50/35]5, 4 inch-wide panel to 
illustrate the mode change which occurred for that panel. 

Recall that each panel was tested two times. The first test was only conducted up to 
the point of failure initiation, so that the location of failure initiation could be determined 
from C-scans, while during the second test, the panels were allowed to fail totally. The 
first three photographs shown for each panel were taken during the first test, while the 
photograph labeled “failure” was taken during the second test and will therefore appear 
slightly different from the other three photographs. 
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Figure B.l. Moire 1 Fringe Photographs for 


Failure 

[±30/90 ] ^ 4 Inch-Wide Panel 
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Before Test 


Buck 1 i ng 



Late Postbuckl i ng 


Fa i lure 



Figure B.2. Moire' Fringe Photographs for [±30/90]^ 5-5 Inch-Wide Panel 
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.ORIGINAL PAGE} IS 
®E POOR QUALITY 




Initial Test Photo 
Unavai lable 


omamM. page „ 

P 00R QUALITY 



Buck 1 i ng 



Late Postbuckl i ng 


Failure 


Figure B.3. Moire 1 Fringe Photographs for [±30/90] ^ 7 Inch-Wide Panel 


1 1 1 


ORIGINAL PAGE IS 

OS SOOR quality 



Late Postbuckl i ng Failure 


Figure B.^. Moire 1 Fringe Photographs for [±50/35] ^ 4 Inch-Wide Panel 
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ORIGINAL PAGE IS' 
OE POOR QUALITY, 



Five Half-Wave Mode Shape Observed Shortly After Buckling 


Six Half-Wave Mode Shape Observed 960 lb Later 



Figure B.5. 

Moire 1 Fringe Photographs Illustrating 
Mode Change Exhibited by [±50/35 ] ^ k Inch-Wide Panel 
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Late Postbuckling Failure 


Figure B.6. Moire 1 Fringe Photographs for [± 50 / 35 ] ^ 5-5 Inch-Wide Panel 
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p AGE is 
-POOR QUALITY 



Before Test 


Buckl i ng 



Late Postbuckl i ng 


Failure 


Figure B.7. Moire' Fringe Photographs for [±50/35] ^ 7 Inch-Wide Panel 
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Appendix C 

C-Scans of Experimental Panels 


This appendix contains C-scans for each of the six tests discussed in Chapter 5. Each 
C-scan is labeled with the location of failure initiation as reported in Figure 5.19, and some 
C-scans which included foreign objects easily confused with damage areas are labeled as 
such. Unfortunately, strain gages and associated wiring axe visible on all C-scans, but 
these objects are easily identifiable. Recall that the white areas indicate damaged areas. 
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gure C.2. C-Scan of [±30/90] 5 5.5 Inch-Wide Panel 
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Of Fai lure 
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Of Failure 


Figure C.3. C-Scan of [±30/90] ^ 7 Inch-Wide Panel 




Figure C.A. C-Scan of [±50/35]^ ** Inch-Wide Panel 
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Figure C.5. C-Scan of [±50/ 35]^ 5*5 Inch-Wide Panel 
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Figure C.6. C-Scan of [±50/35] 5 7 Inch-Wide Panel 
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Appendix D 
STAGSC-1 Runstreams 


This appendix contains the input data corresponding to the STAGSC-1 results pre- 
sented in Chapter 5. The six runstreams differ primarily in the presence or absence of 
constraints (G cards) which force a certain symmetry to occur. Such constraints were 
necessary in order to obtain the postbuckled shape corresponding to the experimentally 
observed shape in cases where the first two eigenvalues were very closely spaced. 
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1 NONLINEAR RUN . . . [+30/-30/90] 5 FLAT PLATE... 4 in 

2 3 1 1 0 0 0 1 $ B1 .'NONLINEAR RUN 

.110001-11 $ B2:l DISPL PARTIAL COMPATIBILITY CARD 

<1010 $ B3:l MATERIAL TABLE; 1 SHELL WALL TYPE 

r, -0.0008 $ B5: FACTOR MULTIPLYING MODE SHAPE USED AS IMPERF. 

6 2000. 500. 8000 . $ Cl .-LOAD FACTOR 

7 0 26000 7-1-1 $ D1:USE ARC LENGTH METHOD 

s 63 11 $ FI: LEVEL OF DISCRETIZATION 

911311101 $ G2: PARTIAL COMPATIBILITY IN U ALONG TOP EDGE 

xo 2 $ G3: SET UP CONSTRAINT TO 

n 1 31 53 1.E7 $ G4: FORCE SYMMETRIC BUCKLE 

12 1 33 73 -1.E7 $ G4: PATTERN 

13 1 0 $ II: MATERIAL TABLE NUMBER 

14 18.5E6 0.03027 0.832E6 1. 1. 1.6E6 0 $ 12: MATERIAL PROPERTIES 

15 1 1 15 0 0 $ Kl:15 LAYERED PLATE 

16 1 5. 11917e-3 -30.0 0 $ K2: LAYER 1 

17 1 5. 11917e-3 30.0 0 $ K2 : LAYER 2 

18 1 5. 11917e-3 -90.0 0 $ K2: LAYER 3 

19 1 5. 11917e-3 -30.0 0 $ K2 : LAYER 4 

20 1 5. 11917e-3 30.0 0 $ K2 : LAYER 5 

21 1 6. 11917e-3 -90.0 0 $ K2 : LAYER 6 

22 1 5. 11917e-3 -30.0 0 $ K2 : LAYER 7 

23 1 5. 11917e-3 30.0 0 $ K2 : LAYER 8 

24 1 5. 11917e-3 -90.0 0 $ K2 : LAYER 9 

25 1 5. 11917e-3 -30.0 0 $ K2: LAYER 10 

26 1 5. 11917e-3 30.0 0 $ K2 : LAYER 11 

27 1 5. 11917e-3 -90.0 0 $ K2 : LAYER 12 

28 1 5. 11917e-3 -30.0 0 $ K2 : LAYER 13 

29 1 5. 11917e-3 30.0 0 $ K2 : LAYER 14 

30 1 5. 11917e-3 -90.0 0 $ K2 : LAYER 15 

31 2 0 $ Ml: SELECT RECTANGULAR GEOMETRY 

32 0. 20. 0. 4. $ M2A: CORNER COORDINATES 

33 1 0 0. 0. 0 0 $ M5 : SHELL WALL RECORD 

34 411 8 6 0 0 $ Nl: 4-NOD ED QUAD ELEMENTS 

35 .375 1.625 6.39 1.61 1.61 6.39 1.625 .375 $ N2: SEGMENTS IN X DIRECTION 

36 1 6 18 6 6 18 6 1 $ N3: ELEMENTS PER SEGMENT 

37 .25 .625 2.25 .625 .25 $ N5:3 SEGMENTS IN Y DIRECTION 

38 1 2 4 2 1 $ N6:# ELEMENTS PER SEGMENT 

39 0 3 2 3 0 $ P1:DEFND BELOW/UNCONSTR/CLAMPED/UNCONSTR 

40 100 000 $ P2:U FREE ON TOP OF PLATE 

41 1 0 0 $ Ql:l LOAD SYSTEM 

42 1 5 0 $ Q2:l LOAD SYSTEM, 5 LOAD CARDS 

43 1.11130 $ q3: APPLIED LOAD ON TOP OF PLATE 

44 0. -1 3 2 0 $ Q3:W=0 ON BOTTOM EDGE OF TOP CLAMPED SUPPORT 

45 0. -1 3 0 10 $ Q3:W=0 ON INSIDE EDGE OF RIGHT SIMPLE SUPPORT 

46 0. -1 3 62 0 $ Q3:W=0 ON TOP EDGE OF BOTTOM CLAMPED SUPPORT 

47 0. -1 3 0 2 $ Q3:W =0 ON INSIDE EDGE OF LEFT SIMPLE SUPPORT 

48 1 0 0 $ Rl: PRINT DISPLACEMENTS 


Figure D.l. STAGS Input Data for [±30/90]s 4 Inch-Wide Panel 
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l HOHLIHEAR RUH . . . [+30/-30/90] 5 FLAT PLATE... 5. 5 in 
23110001 $ Bl: HOHLIHEAR RUM 

310001-11 $ B2:l DISPL PARTIAL COMPATIBILITY CARD 

41010 $ B3:l MATERIAL TABLE; 1 SHELL WALL TYPE 

5 -0.0008 $ B5: FACTOR MULTIPLYIIG MODE SHAPE USED AS IMPERF. 

6 2000. 500. 8000. $ Cl: LOAD FACTOR 

7 0 25000 10 -1 -1 $ D1:USE ARC LEHGTH METHOD 

8 45 13 $ FI : LEVEL OF DISCRETIZATIOI 

911511101 $ G2: PARTIAL COMPATIBILITY II U ALOIG TOP EDGE 

10 2 $ G3: SET UP COISTRAIIT TO FORCE 

u 1 19 63 1.E7 $ G4: AHTISYMMETRIC BUCKLE 

12 1 27 83 1.E7 $ G4: PATTER! 

13 1 0 $ II: MATERIAL TABLE HUMBER 

14 18.5E6 0.03027 0.832E6 1. 1. 1.6E6 0 $ 12: MATERIAL PROPERTIES 

15 1 1 15 0 0 $ Kl:15 LAYERED PLATE 

16 1 5. 16583E-3 -30.0 0 $ K2: LAYER 1 

17 1 5 . 16583E-3 30.0 0 $ K2: LAYER 2 

18 1 5. 16583E-3 -90.0 0 $ K2 : LAYER 3 

19 1 5.16583E-3 -30.0 0 $ K2: LAYER 4 

20 1 5. 16583E-3 30.0 0 $ K2: LAYER 5 

21 1 5. 16583E-3 -90.0 0 $ K2: LAYER 6 

22 1 5. 16583E-3 -30.0 0 $ K2: LAYER 7 

23 1 5. 16583E-3 30.0 0 $ K2 : LAYER 8 

24 1 5. 16583E-3 -90.0 0 $ K2: LAYER 9 

25 1 5. 16583E-3 -30.0 0 $ K2: LAYER 10 

26 1 5. 16583E-3 30.0 0 $ K2 : LAYER 11 

27 1 5. 16583E-3 -90.0 0 $ K2: LAYER 12 

28 1 5. 16583E-3 -30.0 0 $ K2: LAYER 13 

29 1 5. 16583E-3 30.0 0 $ K2: LAYER 14 

30 1 5. 16583E-3 -90.0 0 $ K2: LAYER 15 

31 2 0 $ Ml: SELECT RECTAHGULAR GEOMETRY 

32 0. 20. 0. 5.5 $ M2A : CORKER COORDIHATES 

33 1 0 0. 0. 0 0 $ M5 : SHELL WALL RECORD 

34 41 18500 $ HI :4-HDDED QUAD ELEMEHTS 

35 .376 2.805 4.41 2.41 2.41 4.41 2.806 .375 $ 12 : SEGMEHTS IH X DIRECTIOH 

36 1 696 696 1 * 13: ELEMEHTS PER SEGMEHT 

37 .25 2. 1. 2. .25 $ 15: SEGMEHTS IH Y DIRECTIOH 

38 1424 1 $ H6: ELEMEHTS PER SEGMEHT 

39 0 3 2 3 0 $ PI: DEFID BELOW/UHCOHSTR/ CLAMPED/UHCOHSTR 

40 100 000 $ P2:U FREE OH TOP OF PLATE 

41 1 0 0 $ Ql:l LOAD SYSTEM 

42 1 5 0 $ Q2:l LOAD SYSTEM, 5 LOAD CARDS 

43 1.11150 $ Q3: APPLIED LOAD OH TOP OF PLATE 

44 0. -1 3 2 0 $ Q3:W=0 OH BOTTOM EDGE OF TOP CLAMPED SUPPORT 

45 0. -1 3 0 12 $ Q3:W=0 OH IISIDE EDGE OF RIGHT SIMPLE SUPPORT 

46 0. -1 3 44 0 $ Q3:W=0 OH TOP EDGE OF BOTTOM CLAMPED SUPPORT 

47 0. -1 3 0 2 $ Q3:W=0 OH IISIDE EDGE OF LEFT SIMPLE SUPPORT 

48 1 $ R1 : PRIHT DISPLACEMEHTS OILY 


Figure D.2. STAGS Input Data for [±30/90]s 5.5 Inch-Wide Panel 
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i NONLINEAR RUN. . . [+30/-30/90] S FLAT PLATE... 7. IN 
23110001 $ Bl: NONLINEAR RUN 

31000101 $ B2:l DISPL PARTIAL COMPAT 

41010 $ B3:l MATERIAL TABLE; 1 SHELL WALL TYPE 

5 0.0008 $ B5: FACTOR MULTIPLYING MODE SHAPE USED AS IMPERF. 

6 2000. 500. 8000. $ Cl: LOAD FACTOR 

7 0 99999 10 -1 -1 $ DlrUSE ARC LENGTH METHOD 

8 39 13 $ FI : LEVEL OF DISCRETIZATION 

911511101 $ G2: PARTIAL COMPAT IN U ALONG TOP EDGE 

10 1 0 $ II: MATERIAL TABLE NUMBER 

11 18.5E6 0.03027 0.832E6 1. 1. 1.6E6 0 $ 12: MATERIAL PROPERTIES 

12 1 1 15 0 0 $ Kl:15 LAYERED PLATE 

13 1 5.25083E-3 -30.0 0 $ K2: LAYER 1 

14 1 5 . 25083E-3 30.0 0 $ K2: LAYER 2 

15 1 5.25083E-3 -90.0 0 $ K2: LAYER 3 

16 1 5.25083E-3 -30.0 0 $ K2: LAYER 4 

17 1 5.25083E-3 30.0 0 $ K2: LAYER 5 

18 1 5.25083E-3 -90.0 0 $ K2: LAYER 6 

19 1 5.25083E-3 -30.0 0 $ K2: LAYER 7 

20 1 5.25083E-3 30.0 0 $ K2: LAYER 8 

21 1 6.25083E-3 -90.0 0 $ K2: LAYER 9 

22 1 5.25083E-3 -30.0 0 $ K2: LAYER 10 

23 1 6.25083E-3 30.0 0 $ K2: LAYER 11 

24 1 5.25083E-3 -90.0 0 $ K2: LAYER 12 

25 1 5.25083E-3 -30.0 0 $ K2: LAYER 13 

26 1 5 . 25083E-3 30.0 0 $ K2: LAYER 14 

27 1 5.25083E-3 -90.0 0 $ K2: LAYER 16 

28 2 0 $ Ml: SELECT RECTANGULAR GEOMETRY 

29 0. 20. 0. 7.0 $ M2A: CORNER COORDINATES 

30 1 0 0. 0. 0 0 $ M5: SHELL WALL RECORD 

314118800 $ N1 : 4-NOD ED QUAD ELEMENTS 

32 .375 4.125 2.29 3.21 3.21 2.29 4.125 .375 $ N2 : SEGMENTS IN X DIRECTION 

33 1 84664 8 1$ N3: ELEMENTS PER SEGMENT 

34 .25 .5 1.42 1.33 1.33 1.42 .5 .25 $ N5 : SEGMENTS IN Y DIRECTION 

35 1122 2211$ N6: ELEMENTS PER SEGMENT 

36 0 3 2 3 0 $ PI :DEFND BELOW/UNCONSTR/CLAMPED/UNCONSTR 

37 100 000 $ P2:U FREE ON TOP OF PLATE 

38 1 0 0 $ Ql:l LOAD SYSTEM 

39 1 5 0 $ q2:l LOAD SYSTEM, 5 LOAD CARDS 

40 1.11150 $ Q3: APPLIED LOAD ON TOP OF PLATE 

41 0. -1 3 2 0 $ Q3:W=0 ON BOTTOM EDGE OF TOP CLAMPED SUPPORT 

42 0. -1 3 0 12 $ Q3:W=0 ON INSIDE EDGE OF RIGHT SIMPLE SUPPORT 

43 0. -1 3 38 0 $ Q3: W=0 ON TOP EDGE OF BOTTOM CLAMPED SUPPORT 

44 0. -1 3 0 2 $ Q3:W=0 ON INSIDE EDGE OF LEFT SIMPLE SUPPORT 

45 1 $ Rl: PRINT DISPLACEMENTS ONLY 


Figure D.3. STAGS Input Data for [±30/90] 5 7 Inch-Wide Panel 
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l NONLINEAR RUN. . . [+50/-50/35]5 FLAT PLATE... 4 in 
23110001 $ B1 : NONLINEAR RUN 

310001-12 $ B2:l DISPL PARTIAL COMPAT 

41010 $ B3:l MATERIAL TABLE; 1 SHELL WALL TYPE 

5 0. 0.008 $ B5: FACTOR MULTIPLYING (2ND) MODE SHAPE USED A IMPERF. 

6 2000. 500. 8000. $ Cl: LOAD FACTOR 

7 0 25000 20 -1 -1 $ D1:USE ARC LENGTH METHOD 

8 63 11 $ FI: LEVEL OF DISCRETIZATION 

911311101 $ G2: PARTIAL COMPAT IN U ALONG TOP EDGE 

10 2 $ G3: SET UP CONSTRAINT TO 

11 1 31 53 1.E7 $ G4: FORCE SYMMETRIC BUCKLE 

12 1 33 73 -1.E7 $ G4: PATTERN 

13 1 0 $ II: MATERIAL TABLE HUMBER 

14 18.5E6 0.03027 0.832E6 1. 1. 1.6E6 0 $ 12: MATERIAL PROPERTIES 

15 1 1 15 0 0 $ Kl:15 LAYERED PLATE 

16 1 5. 15417E-3 -60.0 0 $ K2 : LAYER 1 

17 1 5. 15417E-3 50.0 0 $ K2: LAYER 2 

18 1 5. 15417E-3 -35.0 0 $ K2: LAYER 3 

19 1 5. 15417E-3 -60.0 0 $ K2: LAYER 4 

20 1 6. 15417E-3 50.0 0 $ K2: LAYER 5 

21 1 5. 15417E-3 -35.0 0 $ K2: LAYER 6 

22 1 5. 15417E-3 -50.0 0 $ K2: LAYER 7 

23 1 5. 15417E-3 50.0 0 $ K2: LAYER 8 

24 1 5. 15417E-3 -35.0 0 $ K2: LAYER 9 

25 1 5. 15417E-3 -50.0 0 $ K2: LAYER 10 

26 1 5. 15417E-3 50.0 0 $ K2: LAYER 11 

27 1 5. 15417E-3 -35.0 0 $ K2: LAYER 12 

28 1 5. 15417E-3 -60.0 0 $ K2: LAYER 13 

29 1 5. 15417E-3 50.0 0 $ K2: LAYER 14 

30 1 5. 15417E-3 -35.0 0 $ K2: LAYER 15 

31 2 0 $ Ml: SELECT RECTANGULAR GEOMETRY 

32 0. 20. 0. 4. $ H2A: CORNER COORDINATES 

33 1 0 0. 0. 0 0 $ M5: SHELL WALL RECORD 

344118500 $ N1.-4-N0DED qUAD ELEMENTS 

35 .376 1.625 6.39 1.61 1.61 6.39 1.625 .375 $ N2 : SEGMENTS IN X DIRECTION 

36 1 6 18 6 6 18 6 1 $ N3: ELEMENTS PER SEGMENT 

37 .25 .625 2.25 .625 .25 $ N5: SEGMENTS II Y DIRECTION 

38 1 2 4 2 1 $ N6: ELEMENTS PER SEGMENT 

39 0 3 2 3 0 $ P1:DEFHD BELOW/UNCONSTR/CLAMPED/UNCONSTR 

40 100 000 $ P2:U FREE ON TOP OF PLATE 

41 1 0 0 $ qi:l LOAD SYSTEM 

42 1 5 0 $ Q2:l LOAD SYSTEM, 5 LOAD CARDS 

43 1.11130 $ Q3: APPLIED LOAD ON TOP CENTER OF PLATE 

44 0. -1 3 2 0 $ q3:W=0 ON BOTTOM EDGE OF TOP CLAMPED SUPPORT 

45 0. -1 3 0 10 $ q3:W=0 ON INSIDE EDGE OF RIGHT SIMPLE SUPPORT 

46 0. -1 3 62 0 $ q3:W=0 ON TOP EDGE OF BOTTOM CLAMPED SUPPORT 

47 0. -1 3 0 2 $ q3:W=0 ON INSIDE EDGE OF LEFT SIMPLE SUPPORT 

48 1 0 0 $ Rl: PRINT DISPLACEMENTS 


Figure D.4. STAGS Input Data for [±50/35]s 4 Inch-Wide Panel 
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1 NONLINEAR RUN. . . [+S0/-50/35] 6 FLAT PLATE... 5. 5 IN 
23110011 $ B1 .NONLINEAR RUN 

310001-11 $ B2:l DISPL PARTIAL COMPAT 

41010 $ B3:l MATERIAL TABLE; 1 SHELL WALL TYPE 

5 0.008 $ B5 : FACTOR MULTIPLYING MODE SHAPE USED AS IMPERF. 

6 2000. 500. 8000. $ Cl: LOAD FACTOR 

7 0 25000 10 -1 -1 $ D1:USE ARC LENGTH METHOD 

8 45 13 $ FI : LEVEL OF DISCRETIZATION 

911611101 $ G2: PARTIAL COMPAT IN U ALONG TOP EDGE 

10 2 $ G3: SET UP CONSTRAINT TO 

11 1 19 63 1.E7 $ G4: FORCE ANTISYMMETRIC BUCKLE 

12 1 27 83 1.E7 $ G4: PATTERN 

13 1 0 $ II : MATERIAL TABLE NUMBER 

14 18.5E6 0.03027 0.832E6 1. 1. 1.6E6 0 $ 12: MATERIAL PROPERTIES 

15 1 1 15 0 0 $ Kl:15 LAYERED PLATE 

16 1 5. 14333E-3 -50.0 0 $ K2: LAYER 1 

17 1 5. 14333E-3 50.0 0 $ K2: LAYER 2 

18 1 5. 14333E-3 -35.0 0 $ K2: LAYER 3 

19 1 S.14333E-3 -50.0 0 $ K2: LAYER 4 

20 1 5. 14333E-3 50.0 0 $ K2 : LAYER 5 

21 1 5 . 14333E-3 -35.0 0 $ K2: LAYER 6 

22 1 5. 14333E-3 -50.0 0 $ K2: LAYER 7 

23 1 5. 14333E-3 50.0 0 $ K2: LAYER 8 

24 1 5. 14333E-3 -35.0 0 $ K2: LAYER 9 

25 1 5. 14333E-3 -50.0 0 $ K2: LAYER 10 

26 1 5 . 14333E-3 50.0 0 $ K2: LAYER 11 

27 1 5. 14333E-3 -35.0 0 $ K2: LAYER 12 

28 1 5. 14333E-3 -50.0 0 $ K2: LAYER 13 

29 1 5. 14333E-3 50.0 0 $ K2: LAYER 14 

30 1 5. 14333E-3 -35.0 0 $ K2: LAYER 15 

31 2 0 $ Ml: SELECT RECTANGULAR GEOMETRY 

32 0. 20. 0. 5.5 $ M2A: CORNER COORDINATES 

33 1 0 0. 0. 0 0 $ M5 : SHELL WALL RECORD 

344118500 $ Nl: 4-NOD ED QUAD ELEMENTS 

35 .375 2.805 4.41 2.41 2.41 4.41 2.805 .375 $ N2 : SEGMENTS IN X DIRECTION 

36 1 696 696 1 $ N3: ELEMENTS PER SEGMENT 

37 .25 2. 1. 2. .25 $ N5: SEGMENTS IN Y DIRECTION 

38 1424 1 $ N6: ELEMENTS PER SEGMENT 

39 0 3 2 3 0 $ P1:DEFND BELOW/UHCONSTR/ CLAMPED/UNCONSTR 

40 100 000 $ P2:U FREE ON TOP OF PLATE 

41 1 0 0 $ qi:l LOAD SYSTEM 

42 1 5 0 $ Q2:l LOAD SYSTEM, 5 LOAD CARDS 

43 1.11150 $ q3: APPLIED LOAD ON TOP CENTER OF PLATE 

44 0. -1 3 2 0 $ q3:W=0 ON BOTTOM EDGE OF TOP CUMPED SUPPORT 

45 0. -1 3 0 12 $ Q3:W=0 ON INSIDE EDGE OF RIGHT SIMPLE SUPPORT 

46 0. -1 3 44 0 $ Q3:W=0 ON TOP EDGE OF BOTTOM CUMPED SUPPORT 

47 0. -1 3 0 2 $ Q3:W=0 ON INSIDE EDGE OF LEFT SIMPLE SUPPORT 

48 1 $ Rl: PRINT DISPLACEMENTS ONLY 


Figure D.5. STAGS Input Data for [±50/35] 5 5.5 Inch- Wide Panel 
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l NONLINEAR RUN. . . [+50/-50/35]5 FLAT PLATE... 7. IN 
23110001 $ B1 : NONLINEAR RUN 

31000101 $ B2:l DISPL PARTIAL COMPAT 

41010 $ B3:l MATERIAL TABLE; 1 SHELL WALL TYPE 

5 0.0008 $ BS: FACTOR MULTIPLYING MODE SHAPE USED AS IMPERF. 

6 2000. 500. 8000. $ Cl: LOAD FACTOR 

7 0 25000 10 -1 -1 $ D1:USE ARC LENGTH METHOD 

8 39 13 $ FI: LEVEL OF DISCRETIZATION 

911511101 $ G2: PARTIAL COMPAT IN U ALONG TOP EDGE 

10 1 0 $ I 1 : MATERIAL TABLE NUMBER 

11 18.5E6 0.03027 0.832E6 1. 1. 1.6E6 0 $ 12: MATERIAL PROPERTIES 

12 1 1 15 0 0 $ Kl:15 LAYERED PLATE 

13 1 5.20083E-3 -50.0 0 $ K2: LAYER 1 

14 1 5.20083E-3 60.0 0 $ K2: LAYER 2 

15 1 5.20083E-3 -35.0 0 $ K2: LAYER 3 

16 1 5 . 20083E-3 -50.0 0 $ K2: LAYER 4 

17 1 5.20083E-3 50.0 0 $ K2: LAYER 6 

18 1 5.20083E-3 -35.0 0 $ K2: LAYER 6 

19 1 5.20083E-3 -50.0 0 $ K2: LAYER 7 

20 1 5.20083E-3 50.0 0 $ K2: LAYER 8 

21 1 5.20083E-3 -35.0 0 $ K2: LAYER 9 

22 1 5.20083E-3 -50.0 0 $ K2: LAYER 10 

23 1 5. 20083E-3 50.0 0 $ K2: LAYER 11 

24 1 5. 20083E-3 -35.0 0 $ K2: LAYER 12 

25 1 5. 20083E-3 -50.0 0 $ K2 : LAYER 13 

26 1 5.20083E-3 50.0 0 $ K2 : LAYER 14 

27 1 5. 20083E-3 -35.0 0 $ K2: LAYER 15 

28 2 0 $ Ml: SELECT RECTANGULAR GEOMETRY 

29 0. 20. 0. 7.0 $ M2A: CORNER COORDINATES 

30 1 0 0. 0. 0 0 $ M5 : SHELL WALL RECORD 

314118800 $ N1:4-N0DED QUAD ELEMENTS 

32 .375 4.125 2.29 3.21 3.21 2.29 4.125 .375 $ N2 : SEGMENTS IN X DIRECTION 

33 1 84664 8 1$ N3: ELEMENTS PER SEGMENT 

34 .25 .5 1.42 1.33 1.33 1.42 .5 .25 $ N5: SEGMENTS IN Y DIRECTION 

35 1122 2211$ N6: ELEMENTS PER SEGMENT 

36 0 3 2 3 0 $ P1:DEFND BELOW/UNCONSTR/CLAMPED/UNCONSTR 

37 100 000 $ P2:U FREE ON TOP OF PLATE 

38 1 0 0 $ qi:l LOAD SYSTEM 

39 1 5 0 $ Q2 : 1 LOAD SYSTEM, 5 LOAD CARDS 

40 1. 1 1 1 5 0 $ Q3: APPLIED LOAD ON TOP CENTER OF PLATE 

41 0. -1 3 2 0 $ Q3:W=0 ON BOTTOM EDGE OF TOP CLAMPED SUPPORT 

42 0. -1 3 0 12 $ Q3:W=0 ON INSIDE EDGE OF RIGHT SIMPLE SUPPORT 

43 0. -1 3 38 0 $ q3:W=0 ON TOP EDGE OF BOTTOM CLAMPED SUPPORT 

44 0. -1 3 0 2 $ Q3:W=0 ON INSIDE EDGE OF LEFT SIMPLE SUPPORT 

45 1 $ Rl: PRINT DISPLACEMENTS ONLY 


Figure D.6. STAGS Input Data for [±50/35]s 7 Inch-Wide Panel 
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Deformed Configuration 
Undeformed Configuration 



Figure 2 . 2 . Effect of Shear-Extens ional Coupling 
On Panel Subjected to Uniform Axial 
Load i ng 


Note Skewing of Nodal Lines 



Figure 2.3. Effect of Shear-Extensional Coupling on Postbuckled Shape 
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Load 


1 2 



Location of i Displacement 
Transducers in Reflectionally 
Symmetric Positions 


Typical Results From Such 
Measurements 


Figure 2 , 4 . Method for Observing Twisting 
Behavior During an Experiment 
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Transverse Loading Problems: 



Axial Loading Problems: 



Figure 4.1. Schematics for Linear Problems 
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Half the panel is modeled: 


Node group 1 


u 2 = = ^2 = 


Node group 2 


0 


t Center of Panel 


b 


For axial loading 
problems, 

Uj = constant 
along this edge. 


a 

2 


1. Node group 1 is related to node group 2 by the symmetry 
relations of Chapter 2 and is therefore not modeled. 

2. Left half of panel is related to right half by the 
symmetry relations of Chapter 2 and is not explicitly 
modeled . 

3. Boundary conditions are imposed as described in the 
previous figure for the three exterior edges, except 
as noted in the above figure. 


Figure b.2. Schematic of Mixed Formulation Finite Element Models 



max 


w max = 4 . 80X1 0 -2 in. 




Normalized Distance Along Long Panel Centerline 


LEGEND 

o Mixed Formulation Finite Element Solution 

Rayleigh Ritz Solution, 3 Terms 

Rayleigh Ritz Solution, 9 Terms 

Rayleigh Ritz Solution, 19 Terms 

Rayleigh Ritz Solution, 33 Terms 

Rayleigh Ritz Solution, 51 Terms 

• Only w and <f> j are shown since u and v are zero 
at all points, and <f> 2 is zero along the long centerline. 


Figure 4.3. Linear Results for Simply Supported 
Transversely Loaded [04/904] s Panel 
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Figure 4.4. Linear Results for Simply Supported 
Transversely Loaded [±50/35]5 Panel 
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Normalized Distance Along Panel Centerlines 


LEGEND 

o Mixed Formulation Finite Element Solution 

Rayleigh Ritz Solution, 9 Terms 

Rayleigh Ritz Solution, 19 Terms 

Rayleigh Ritz Solution, 33 Terms 

Rayleigh Ritz Solution, 51 Terms 

• u is plotted down the long panel centerline, 
and v is plotted down the short centerline. 

• Only u and v are shown since all other 
displacements are zero at all points. 


Figure 4.5. Linear Results for Clamped/Simply Supported 
Axially Loaded [0 4 /90 4 ] s Panel 
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Figure 4.6. Linear Results for Clamped/Simply Supported 
Axially Loaded [±50/35]5 Panel 
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ui Displacement 



ui Displacement 


Figure 4.7. Contour Plots of Linear Solution For Axially Loaded 
Clamped/Simply Supported [±50/35]s Panel 


139 







LEGEND 

a Rayleigh Ritz Solution, 19 Terms 

b Rayleigh Ritz Solution, 33 Terms 

c Mixed Formulation Finite Element Solution 

d STAGS Solution 


Figure 4.8. End Shortening Behavior, [04/90 4 ] s Panel 
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LEGEND 

a- 

Rayleigh Ritz Solution, 19 Terms 

b- 

Rayleigh Ritz Solution, 33 Terms 

c- 

Mixed Formulation Finite Element Solution 

d- 

STAGS Solution 


Figure 4 . 9 . Center Transverse Displacement Behavior, [04/904]$ Panel 
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Rayleigh-Ritz Results at P = 4407 lb 


Figure 4.10. Normalized Contour Plots of Transverse 
Displacement, [O 4 / 9 O 4 ] Panel 
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End Shortening, In. 


LEGEND 

a Rayleigh Ritz Solution, 38 Terms 

b - Rayleigh Ritz Solution, 66 Terms 

c Mixed Formulation Finite Element Solution 

d — - — STAGS Solution 
e — Experiment 


Figure 4 . 11 . End Shortening Behavior, [±30/90)5 Panel 
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LEGEND 

a Rayleigh Ritz Solution, 38 Terms 

b Rayleigh Ritz Solution, 66 Terms 

c Mixed Formulation Finite Element Solution 

d STAGS Solution 

e Experiment 

Figure 4.12. Center Transverse Displacement Behavior, [=b30/90]s Panel 
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Rayleigh-Ritz Results at P = 7137 lb 


Figure 4.13. Normalized Contour Plots of Transverse 
Displacement, [±30/90]s Panel 
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Load, P/P 


4 . 


0.06 0.08 0.10 
End Shortening, In. 


LEGEND 

a Rayleigh Ritz Solution, 38 Terms 

b - Rayleigh Ritz Solution, 66 Terms 

c Mixed Formulation Finite Element Solution 

d STAGS Solution 

e Experiment 

Figure 4.14. End Shortening Behavior, [±50/35]s Panel 
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W, In 


LEGEND 

a Rayleigh Ritz Solution, 38 Terms 

b - - • - - Rayleigh Ritz Solution, 66 Terms 

c Mixed Formulation Finite Element Solution 

d STAGS Solution 

e — - — Experiment 

Figure 4.15. Center Transverse Displacement Behavior, [±50/35]s Panel 
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Rayleigh-Ritz Results at P = 6325 lb 


Figure 4.16. Normalized Contour Plots of Transverse 
Displacement, [±50/35]s Panel 
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a. Schematic of Experiment Setup 



b. Photograph of Experiment Setup 
Figure 5-1. Experiment Schematic and Photograph 
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Loaded Edge 


a. Schematic of Clamping Fixtures 


Knife-Edge 

Simple Support Fixtures 



Section B-B 


Loaded Edge 

b. Schematic of Simple Support Fixtures 
Figure 5 . 2 . Support Fixture Schematic 
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Load 
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Coupl ing 


Load 
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Degree of Bending/ 
Extensional Coupling 


Figure 5.3. 


Tension 



Deformed 

Configuration 




Bending as Indicated By 
Back-to-Back Strain Gage Data 
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Figure 5*^- Sample C-Scan 
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Displacement 
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X 
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////// 

Y / / / / / 


Straight 

Edge 


To ensure that load is introduced parallel to the 
plane of the panel : 

- Isotropic Panel: Points (1) and (3) must exhibit 
identical strain. (These points are reflectionally 
symmetric.) 

- Anisotropic Panel: Points (1) and (2) must exhibit 
identical strain. (These points are inversionally 
symmetric.) 

Figure 5.5, Comparison of Experimental 

Setup Technique for Isotropic 
And Anisotropic Panels 
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Figure 5.6. Undeformed Geometry Plot of STAGSC-1 Finite 
Element Model for 5.5" [±50/35]5 Panel 
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Load, lb X 1 0 3 



Figure 5.7. Typical Panel End Shortening Plot 
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a. [±30/90], 7 In. f. [±.50/35], 7 ?n. 

End Shortening, in. 


LEGEND 

- STAGS Nonlinear Analysis 

Experiment 

• Horizontal lines indicate estimated buckling loads. 

Figure 5.8. End Shortening Results 
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Load, lb X 10 3 



LEGEND 

STAGS Nonlinear Analysis 

Experiment 

• Calculation is made at loads of 2P C r for both curves. 

• Data shown is for [±50/35] 5 5.5 inch-wide panel. 

• Horizontal lines indicate esti mated buckling loads. 

Figure 5.9. Estimation of Stiffness Degradation 
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Near-Maximum Transverse Displacement, in. 


LEGEND 

STAGS Nonlinear Analysis 

Experiment 

• Horizontal lines indicate estimated buckling loads. 

Figure 5.10. Near-Maximum Transverse Displacement Results 
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0*00 0.0S 0.10 0.15 0.20 0.25 0.30 

a. [±30/90] s 4 In. 


0.00 0.05 0.10 0.15 0.20 0.25 0.30 

b. [±50/35] 3 4 In. 
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r 

0 1 1 ■ ■ 1 1 1 ■ ■ 1 1 ■ i ■ i ■ ■ ■ ■ i ■ ■ 1 1 i-i 
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c. [±30/90], 5.5 In. 


' ■ ' i 
0.30 



d. [±50/35], 5.5 In. 




«. [±30/90] s 7 In. 


f. [±50/35], 7 In. 


Absolute Value of Transverse Displacement, in. 


LEGEND 


ft 

ft 

• Displacements are plotted at two 
reflectionally symmetric points as shown. 



• Horizontal lines indicate estimated buckling loads. 


Figure 5.11. Experimental Twisting Results 
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43 



4. 

Load, 

2. 



0.00 0.05 0.10 0.15 0.20 0.25 0.30 

o. [±30/90] s 5.5 In. 




0.00 0.05 0.10 0.15 0.20 0.25 0.30 

d. [±50/35]. 5.5 In. 



e. [±30/90] # 7 In. f. [±50/35] s 7 In. 

Absolute Value of Transverse Displacement, in. 

LEGEND 


• Displacements are plotted at two 
reflectionally symmetric points as shown. 


• Horizontal lines indicate estimated buckling loads. 


Figure 5.12. STAGSC-1 Twisting Results 


160 








£ 0 i 


0.00 0.05 0.10 0.15 0.20 0.25 0.30 

a. [±30/90], 4 In. 


0.00 0.05 0.10 0.15 0.20 0.25 0.30 


b. [±50/35], 4 In. 



0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.00 0.05 0.10 0.15 0.20 0.25 0.30 


«. [±30/90], 7 In. f. [±50/35], 7 In. 

Absolute Value of Transverse Displacement, in. 



Figure 5.13. Experimental Symmetry Results 
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c. [±30/90], 5.5 In. d. [±50/35], 5.5 In. 




a. [±30/90], 7 In. 

Absolute Value of Transverse Displacement, in. 


LEGEND 



* 

a Displacements are plotted at two 


* 


inversionally symmetric points as shown. 

a 

Horizontal lines indicate estimated buckling loads. 


Figure 5.14. STAGSC-1 Symmetry Results 


162 



a. [±30/90] s 4 In. 



b. [±50/35] 3 4 In. 


to 
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_J 

- 0.010 - 0.003 0.000 0.003 0.010 

c. [±30/90] s 5.5 in. 




d. [±50/35],, 5.5 in. 




f. [±50/35] s 7 In. 

Representative Surface Strains, in/in 


LEGEND 

• For 4 and 7 inch-wide panels, center panel strains are shown. 

• For 5 inch-wide panels, strains near top of panel in the 
center of the top half-wave are shown. 

• Horizontal lines indicate estimated buckling loads. 

Figure 5.15. Representative Experimental Surface Strain Results 
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Back-to-Back Strains Corresponding to Points Noted Below 



Distance Down Long Centerline of Panel, in. 

Figure 5.16. Progression of Transverse Displacement Pattern, [±50/35] 5 7 Inch-Wide Panel 






Back-to-Back Strains Corresponding to Points Noted Below 



Figure 5.17. Progression of Transverse Displacement Pattern, [±30/90]5 7 Inch-Wide Panel 





Transverse Displacement, in. 




e. [±30/90], 7" (P cr = 2500 lb) f. [±50/35], 7 " (P or = 2670 lb) 

Distance Down Long Centerline of Panel, in. 


LEGEND 

• Symbols indicate the inflection points on each curve, 
as determined using cubic spline interpolation. 

Figure 5.18. Progression of Transverse Displacement Pattern, All Panels 
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Figure 5.19. Failure Locations 
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• Experimental surface strain values 
o Calculated average strain values 
Reference lines to aid in observing symmetry 

Figure 5.20. Strains Across 5.5 Inch-Wide Panels at Load Near Failure 
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